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1.  INTRODUCTION 


The  description  of  the  earth's  gravity  field,  or  in  other  wordB,  its 
approximation  with  mathematical  models,  is  an  old  interest  of  physical  geodesy 
and  more  recently  it  is  a  common  goal  with  geophysics  in  general  and  with 
satellite  orbit  determinations  in  particular.  Before  the  development  of  satellite 
techniques  the  observation  techniques  were  restricted  to  terrestrial 
gravimetry;  adequate  gravity  coverage  existed  only  over  industrially  well- 
developed  areas.  The  gravity  profiles  over  the  oceans  were  inadequate  in 
coverage,  inefficient,  and  expensive,  in  addition  to  doubtful  accuracy  of  the 
results.  Under  these  circumstances,  the  development  of  a  meaningful  gravity 
model  was  impossible. 

The  satellite  tracking  technology  opened  the  potential  of  geopotential 
modeling  by  addition  of  new  types  of  efficient  data  collection  techniques  with 
global  coverage  capabilities.  With  the  improvements  in  the  quality,  quantity, 
and  distribution  of  satellite  tracking  data  (Smith  1982);  as  well  as  the  increase 
in  coverage,  refinements  in  resolution,  and  improvements  in  accuracy  of 
surface  data  (Rapp  1978,  1981,  1983;  Rowland  1981);  the  geopotential  models 
showed  significant  progress.  Satellite  tracking  contributed  to  the  long 
wavelength  part  of  the  spectrum  expressed  by  improved  low  degree 
coefficients;  satellite  altimetry  and  improved  terrestrial  surface  information 
contributed  to  the  medium  lengths  of  the  spectrum.  New  developments  in  data 
analysis  and  combination  techniques  helped  to  produce  improved  geopotential 
models. 

Despite  the  improvements,  large  discrepancies  still  exist  between  the 
various  models,  even  between  those  derived  from  nearly  the  same 
observations.  Some  are  "tailored"  for  specific  satellite  orbits,  others 
represent  better  surface  data  in  specific  areas.  A  general  model  which  is 
equally  adequate  for  orbit  determinations  and  geodetic  use  is  lacking  at  the 
present  time;  the  principal  reason  is  the  lack  of  observation  material  covering 
the  full  spectral  range. 

In  the  preceding  years  programs  were  organized  by  various  institutions 
both  in  the  United  States  and  in  Europe  for  a  systematic  approach  with 
specific  goals  for  the  improvement  of  the  geopotential.  The  Committee  on 
Geodesy  of  the  National  Research  Council,  National  Academy  of  Sciences,  U.S.A., 
in  "Applications  of  a  Dedicated  Gravitational  Satellite  Mission",  Washington  D.C., 


1979  suggested  research  goals  for  the  geodynamic  program  of  NASA.  These 
include  a  gravitational  satellite  mission  (GRAVSAT)  to  provide  data  necessary 
for  the  achievement  of  the  objectives  of  the  program  (Lerch,  1983).  The 
accuracies  required  with  resolutions  at  100  to  3000km  were  2.5  to  lOmgal  for 
gravity  anomalies  and  10cm  accuracy  for  the  oceanic  geoid.  Later  NASA 
developed  a  "Geopotential  Research  Program  Plan"  (NASA  1982)  with  specific 
goals  for  gravity:  *lmgal  accuracy  with  100km  resolution  for  the  gravity 
field;  5cm  accuracy  with  100km  resolution  for  the  geoid  (Murphy,  1983).  This 
program  is  divided  into  three  phases:  a)  interim  field  model  improvements, 
b)  geopotential  research  mission,  and  c)  advanced  mission.  Currently  all  three 
phases  are  progressing  simultaneously  and  they  will  be  discussed  later  in  this 
report. 

The  European  Space  Agency  (ESA)  also  prepared  a  plan  called  Space  Laser 
Low  Orbit  Mission  (SLALOM)  to  obtain  global  data  coverage  for  improved 
geopotential  modeling. 

Current  geopotential  models  contain  the  low  and  medium  frequency  parts 
of  the  spectrum,  because  the  data  types  (satellite  tracking  data  and  1*  ■  1* 
mean  anomalies)  do  not  contain  high  and  very  high  frequency  information. 
Theoretically  each  data  type  contains  the  total  frequency  range,  however  in 
reality  the  measuring  process  acts  as  a  bandpass  filter  limiting  the  range  of 
the  spectrum  (Schwarz,  1984).  Therefore,  the  combination  of  different  types 
of  measurements  is  necessary  to  obtain  a  homogeneous  spectral  resolution. 
Data  types  containing  high  frequency  are:  5’  x  5’  anomalies,  deflections  of  the 
vertical,  inertial  survey  data,  and  airborne  gradiometry.  The  very  high 
frequency  field  can  be  obtained  from  point  gravity  measurements,  inertial 
surveys,  airborne  gradiometry,  height  data,  and  surface  density  information. 

The  available  high  frequency  gravity  data  covering  a  limited  area  are 
usually  used  for  the  determination  of  some  functionals  of  the  anomalous 
potential  in  a  local  area  or  at  some  selected  points.  These  functionals  can  be 
expressed  mathematically  by  integral  formulae  (space  domain)  or  by  spherical 
harmonic  series  (frequency  domain).  Theoretically  the  two  types  of  formulae 
are  equivalent,  but  in  reality  they  are  different  because  of  the  differences  in 
the  types  and  other  characteristics  of  the  gravity  information  used  for  the 
computations.  The  quality  of  the  estimation  of  a  local  or  regional  gravity  field 
depends  on  several  factors  such  as:  a)  the  sensitivity  of  the  gravimetric 
quantity  to  be  estimated  to  the  given  type  of  data  set,  b)  the  density  of  the 


data  set,  and  c)  area  of  the  coverage  and  measurement  accuracy.  Various 
functionals  are  sensitive  to  different  frequency  ranges  represented  by 
different  types  of  data:  for  example,  in  accurate  geoid  determination  the  low 
frequency  data  are  dominant;  in  accurate  determination  of  the  radial 
component  of  the  second  order  gravity  gradient  the  high  and  very  high 
frequency  data  are  required. 

Recently  local  gravity  data  Bets  and  elevation  data  are  available  in 
gridded  form  and  high  degree  geopotential  models  can  be  used  as  reference 
fields.  As  a  result,  the  higher  frequency  features  of  the  gravity  field  can  be 
analyzed  in  local  areas.  The  problem  in  this  case,  as  in  the  case  of 
combination  of  global  data  sets,  is  the  weighting  of  data  sets  in  a  combination 
solution. 

In  the  following  sections  of  this  report  the  approximation  of  the  external 
gravity  potential  of  the  earth  is  discussed  in  detail.  The  description  of  the 
global  gravity  field  by  spherical  harmonic  series  includes  most  of  the  gravity 
models  published  since  the  late  1970’s.  The  observation  programs  and 
techniques  utilized  now  and  in  the  planning  Btage  for  future  use  are 
described.  The  present  techniques  and  potential  improvements  for  adjustment 
and  data  combination  for  gravity  modeling  are  summarized. 

Regarding  the  local  and  regional  gravity  field  estimation  the  following 
topics  are  discussed: 

(a)  the  spectral  properties  of  various  data  types 

(b)  estimation  of  gravimetric  quantities  in  local  and 
regional  areas. 

The  estimation  techniques  covered  by  the  reviews  of  several  recent 
studies  include  the  use  of  modified  integral  formulae  and  of  the  collocation 
method. 

Due  to  the  specific  interest  in  this  topic,  the  prediction  of  gravity 
disturbance  at  high  altitude  is  also  covered. 

Considering  that  the  high  and  very  high  frequency  local  and  regional 
field  improvement  is  expected  from  airborne  gradiometry,  several  accuracy  and 
feasibility  studies  in  this  area  are  also  reviewed. 


2.  GEOPOTENTIAL  MODBLS 


The  gravitational  potential  of  the  earth  can  be  expressed  in  an  earth  fixed 
and  earth  centered  coordinate  system  by  the  well  known  harmonic  series: 


V(r,4,X)  =  2! 


^  M 

cos  n»X  +  S tm  sin  mxj  *  Pgm  (sin  ♦) 


i=2  m=o 


where: 


r,4,X 


are  the  geocentric  coordinates  of  a  point. 


04  is  the  product  of  the  gravitational  constant  (G)  and  the 

aass  of  the  earth  (M) , 

a  is  a  scale  factor,  usually  the  equatorial  radius  of  the 

reference  ellipsoid, 

C|M, Sj|n  are  a  set  of  fully  normalized  harmonic  coefficients. 


are  fully  normalized  Legendre  functions. 


The  coefficients  Cin  and  SjM  are  then  determined  up  to  a  certain  degree  and 
order  from  analysis  of  satellite  orbit  perturbations  with  or  without  combination 
of  surface  gravity  data. 

In  the  following  sections  some  of  the  recent  contributions  to  the  global 
geopotential  model  development  will  be  recapitulated  in  groups  arranged 
according  to  the  originating  institutions. 


2.1  Goddard  Barth  Models  (GEM)  of  the  Goddard  Space  Plight  Center  (NASA), 
Greenbelt,  Maryland,  20771 

2.1.1  GBM  9  and  GEM  10  Models 

GEM9  (Lerch  and  al.  1979)  was  derived  by  using  laser  tracking  data  from 
GEOS-3,  LAGEOS  and  Starlette  satellites;  S-band  measurements  on  LANDSAT  I; 
and  various  tracking  data  on  26  other  satellites  used  in  previous  GEM  models. 

These  represent  about  840,000  satellite  tracking  measurements,  200,000  of 
which  are  laser  ranging  data.  Because  of  the  sensitivity  of  laser  ranging  to 
satellite  perturbations  observation  partials  on  nine  satellites  were  computed 


complete  through  degree  and  order  22  for  the  harmonic  coefficients.  For 
other  satellites,  measurements  are  complete  only  to  degree  and  order  16.  The 
orbital  data  and  the  list  of  satellites  observed  are  given  in  Lerch  and  al. 
(1979).  GEM9  harmonic  coefficients  are  complete  to  order  and  degree  20  with 
selected  higher  degree  terms. 

GEM10  (Lerch  and  al.  1979)  contains  surface  gravity  data  in  the  form  of 
5*  «  5*  mean  anomalies  (Rapp  1977).  Of  the  1654  5*  *  5*  blocks,  1507  are 
computed  from  1*  *  1*  values  within  each  block,  and  174  5*  *  5*  mean  values 
were  interpolated.  These  surface  date  were  combined  with  the  satellite 
information  of  GEM9.  GEM10  is  complete  to  degree  and  order  22  with  selected 
higher  degree  terms  to  degree  and  order  30. 

In  the  development  of  GEM9  and  GEM10,  changes  were  made  in  the 
modeling  technique  used  for  previous  GEM  solutions.  One  of  these  changes 
was  the  use  of  a  "modified  least  squares  solution"  (Lerch  and  al.  1979).  By 
using  Kaula’s  rule  ( aa  -  10“*/la)  for  the  coefficients  as  a  constraint,  the 
extension  of  the  solution  to  order  and  degree  20  was  achieved.  The  use  of 
least  squares  collocation  not  only  permitted  a  larger  field  but  the  application 
of  least  squares  collocation  was  also  required  to  achieve  acceptable  results 
due  to  the  high  correlation  between  several  of  the  higher  degree  and  order 
coefficients  (Lerch  and  al.  1981).  Another  modification  was  the  use  of  less 
weight  for  the  surface  gravity  data  in  GEM  10. 

Geocentric  coordinates  for  146  tracking  stations  were  also  determined  with 
estimated  accuracies  of  l-2m  in  each  direction.  The  equatorial  radius  of  the 
reference  ellipsoid  was  derived  by  three  different  methods,  all  agreeing  within 
lm  to  6378139m;  the  value  for  GM  was  also  determined  from  laser  ranging  as 
GM  =  398,600.64  km3/seca.  The  value  of  the  flattening  1/f  =  298.257  (updated 
by  the  recovery  of  Ja). 

The  GEM9  and  10  gravity  fieldB  were  evaluated  for  orbit  determination 
accuracy  and  compared  with  other  independent  information  for  geoid  heights 
and  gravity  anomalies.  GEOS-3  orbital  accuracies  are  about  lm  in  radial 
components  for  5  day  arc.  The  accuracies  of  the  coefficients  imply  commission 
errors  in  geoid  height  *1.9m  RMS  for  GEM9  and  *1.5m  RMS  for  GEM10. 

The  harmonic  coefficients,  station  coordinates,  results  of  comparisons  and 
evaluations  are  tabulated  or  plotted  in  the  paper:  Lerch  and  al.  (1979). 


GBMIOB  (Larch  and  al.  1981)  ia  a  model  derived  by  the  combination  of 
GBMIO  data  with  about  700  altimetry  passes  of  GEOS-3  globally  distributed 
with  a  2*  spacing.  GBMIOB  solution  was  a  least  squares  adjustment  without 
any  constraints;  it  is  complete  to  order  and  degree  36.  The  geoid  was  also 
computed  from  the  harmonic  coefficients  using  B  run’s  formula. 

GBM10C  is  an  extension  of  GBMIOB  beyond  degree  and  order  36  and  it  ia 
complete  to  degree  and  order  180;  it  has  the  same  harmonic  coefficients  as 
GBMIOB  through  degree  36.  The  extended  field  was  derived  from:  a)  38,000 
1*  *  1*  mean  gravity  anomalies  from  surface  measurements  (Rapp  1979a)  (The 
geophysically  predicted  surface  anomalies  in  the  Rapp  file,  predicted  by  DMA 
Aerospace  Center,  were  not  used);  b)  GEOS-3  altimeter  data  in  the  form  of 
28,000  1*  *  1*  mean  sea  surface  values  (geoid  heights).  Over  3000  passes 
have  been  used  to  produce  the  set  of  mean  surface  values.  In  combining  the 
data,  a  single  value  per  1*  block  was  used  and  altimeter  data  was  preferred 
over  surface  measured  values.  In  total,  about  50,500  1*  ■  1*  blocks  with 
gravity  or  altimeter  observations  were  used  in  the  adjustsient  of  GBM10C. 
Residual  geoid  heights  were  computed  by  subtraction  of  GBMIOB  geoid  from 
the  1*  global  set  of  geoid  heights.  The  computer  program  developed  by  Risos 
(1979)  was  used  for  the  harmonic  analysis. 

The  use  of  altimeter  data  resulted  in  significant  improvements  in  the  GEM 
models.  Extensive  tests  have  been  made  to  determine  the  improvements  in 
terms  of  radial  orbit  and  geoid  accuracies. 

The  radial  accuracy  of  the  modelB  were  evaluated  by  computing  GEOS-3 
orbits  from  the  same  laser  ranging  data  using  different  gravity  siodels.  The 
difference  of  the  altimeter  residuals  at  the  crossover  points  is  a  measure  of 
the  radial  orbital  error.  A  crossover  test  comparing  GBM10  and  GBMIOB  using 
348  altimeter  crossover  points  resulted  in  orbital  RMS  errors  for  GBM10  of 
1.37m  and  for  GBMIOB  of  1.00m. 

To  teat  the  accuracy  of  the  global  geoids  obtained  from  these  models,  the 
RMS  fit  between  the  GEOS-3  altimeter  profiles  and  those  computed  from  the 
models  are  listed  for  several  GEM  models  in  Lerch  and  al.  (1981).  Residual 
RMS  of  fit  for  10  "trench  arcs"  and  for  10  "non-trench  arcs"  are  respectively: 
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Only  a  little  improvement  can  be  seen  for  the  "trench"  passes  until  high 
degree  coefficients  are  computed  (GEM10C).  For  the  "non-trench"  passes,  the 
improvement  is  100%  in  GEM10B  versus  GEM10.  SEASAT  profiles  were  also  used 
for  comparison  with  GEM10B  values.  The  two  geoids  match  within  *lm. 

The  geodetic  parameters  computed  from  GEM10B  are  very  close  to  current 
values  adopted  by  the  LAG:  GM  =  398600.44  *  .02  kmVsec2,  a.  =  6378138  *  lm, 
1/f  =  298.257  *  0.001,  and  g.  =  978031.5  *  0.3  mgal. 

Considering  that  GEM10C  is  an  extension  of  GEM10B  by  harmonic  analysis 
of  the  residual  1*  *  l*  geoid  heights,  the  authors  (Lerch  and  al.  1981) 
compared  these  residual  geoid  heights  with  values  computed  from  the  harmonic 
coefficients  of  the  "residual"  gravity  field.  The  RMS  difference  was  44cm, 
which  is  a  very  good  agreement.  The  power  spectrum  of  the  180  degree  field 
"shares  an  excellent  continuity"  with  the  GEM10B  spectrum  for  terms  beyond 
degree  36. 


2.1.3  GEM  L-2  model 

A  new  Goddard  Earth  Model,  GEM  L-2,  was  developed  from  "satellite  only" 
observations  utilizing  2  1/2  years  of  LAGEOS  satellite  laser  ranging  data  in 
combination  with  tracking  data  on  30  other  satellites  in  GEM9  (Lerch  and  al. 
1982a).  The  shape,  high  density  and  high  altitude  (about  one  earth  radius)  of 
LAGEOS  eliminates,  for  all  practical  purposes,  errors  due  to  uncertainties  in 
modeling  non-gravitational  forces.  This  model,  complete  through  degree  and 
order  20,  contains  more  that  600,000  laser  measurements  with  more  than  half  of 
them  on  LAGBOS.  Due  to  the  high  altitude  of  LAGEOS,  the  long  wavelength 
features  are  substantially  improved;  however,  the  coefficients  beyond  degree  7 
are  dominated  by  GEM9  data.  The  main  purpose  of  this  model  was  to  obtain 
more  precise  station  positions  and  in  turn,  better  distances  between  the 
stations  for  NASA's  crustal  dynamics  program  and  to  improve  the  long 
wavelength  gravity  field.  The  same  modified  least  squares  method  that  was 


employed  for  the  GEM9  model  was  utilized  for  this  solution.  The  method  is 
described  in  Lerch  and  al.  (1977) 

Simultaneously  with  the  harmonic  coefficients  and  station  coordinates,  the 
GM  term  was  computed  as  398600.607  kmVsec*,  using  c  =  299792.5  m/sec  as  the 
speed  of  light.  Polar  motion  and  Al-UTl  variations  were  estimated  from  5-day 
segments  of  LAGEOS  data.  Finally  earth  tides  were  modeled  with  Love 
numbers  h  =  .60,  i2  ~  .075,  and  k  =  .29.  The  solid  earth  phase  angle  was 
obtained  as  2.018*. 

The  improved  accuracies  of  GEM  L-2  coefficients  versus  those  of  the  GEM9 
model  through  degree  and  order  four  are  significant.  They  are  tabulated  in 
(Lerch  and  al.  1982a).  An  independent  verification  of  these  accuracies  has 
been  performed  by  Wagner  (1982)  using  independent  long  term  24  hour 
satellite  longitude  accelerations.  GEM  L-2  satisfies  these  accelerations.  The 
'by  degree’  contribution  of  the  harmonic  coefficients  to  the  total  longitude 
accelerations  show  that:  "while  the  accelerations  are  sensitive  to  harmonics  as 
high  as  degree  6,  the  GEM  L-2  uncertainties  beyond  degrees  3  or  4  are  not 
directly  tested  by  these  longitude  accelerations"  (Lerch  and  al.  1982a). 

The  improved  GEM  L-2  resulted  in  more  accurate  orbits,  and  in  turn,  more 
accurate  station  positions.  Consequently,  the  capabilities  to  improve 
monitoring  tectonic  motions  by  satellite  laser  ranging  were  also  improved.  The 
polar  motion  and  Al-UTl  variation  estimates  contributed  to  the  stability  of 
harmonic  coefficient  solutions  and  improved  positioning  accuracy.  The  polar 
motion  values  have  a  precision  of  about  10  cm. 

The  effects  of  the  GEM9  and  GEM  L-2  models  and  of  the  polar  motions 
from  B1H  and  GEM  L-2  on  station  positioning  were  compared.  For  the 
comparison  the  entire  LAGEOS  data  set  (1979  through  1981)  has  been  divided 
into  two  independent  sets  by  splitting  each  month’s  data  into  two  15-day 
segments.  This  way  the  two  solutions  contain  independent  data,  span  the 
same  time  interval,  and  average  out  any  resulting  plate  motion  over  these  base 
lines  (distances).  Comparison  of  28  distances  between  eight  stations  were 
made  using  GEM9  and  GEM  L-2  station  coordinates  and  with  GEM  L-2  and  B1H 
polar  motion.  The  RMS  error  for  the  28  lines  from  GEM9  was  7.2  cm;  the  same 
from  GEM  L-2  model  and  polar  motion  was  only  1.8  cm. 

In  summary,  the  laser  ranging  on  LAGEOS  and  the  GEM  L-2  model  resulted 
in:  a)  an  improved  method  for  tracking  station  positioning,  b)  a  great 

improvement  of  the  harmonic  coefficients  of  the  gravity  field  through  degree 


and  order  4,  and  c)  an  improvement  of  LAGEOS  orbit  errors  from  lm  with 
GBM9  to  about  30  cm  with  GEM  L-2  (F.J.  Lerch  and  al.  1982a,  1984). 


2.2  The  Ohio  State  University  (Rapp)  Models 
2.2.1  The  Rapp  1978  Model 

In  this  combination  solution  the  following  satellite,  terrestrial  gravity,  and 
altimeter  data  were  used  (Rapp  1978): 

—  The  potential  coefficients  of  GEM9  (Lerch  and  al.  1979)  derived  from 

satellite  tracking  data  up  to  degree  and  order  20  were  used  in  two 
adjustments:  the  first  adjusted  the  coefficients  up  to  degree  8  and  the 

second  to  degree  12. 

—  The  second  type  of  data  was  a  Bet  of  1*  *  1*  mean  anomalies  consisting  of 
39,405  land  anostalies  and  29,478  1*  *  1*  mean  values  derived  from  GEOS-3 
altimeter  data  (Rapp  1978a).  After  editing  and  replacing  oceanic  estimated 
values  with  altimeter  anomalies,  where  available,  a  total  of  50,650  1*  ■  1*  block 
values  were  retained  in  the  merged  set.  For  the  miBsing  14150  anomalies,  zero 
values  were  used  with  a  standard  deviation  of  *30  mgal. 

The  combination  method  essentially  is  a  weighted  adjustment  of  the 
harsionic  coefficients  computed  from  the  mean  anomalies  and  those  obtained 
from  the  satellite  model.  The  results  are  a  set  of  consistent  harmonic 
coefficients  and  mean  anomalies.  Details  of  the  method  are  given  in  Kaula 
(1966),  Rapp  (1968),  Snowden  and  Rapp  (1968),  and  Rapp  (1978). 

Each  solution  resulted  in  the  adjusted  harmonic  coefficients  (to  degree  8 
and  12  respectively)  and  a  set  of  64800  1*  *  1*  adjusted  mean  anomalies. 
Bach  set  of  the  mean  anomalies  were  developed  into  harmonic  coefficients  to 
degree  30  and  compared  to  the  coefficients  of  the  GEM9  solution.  The 
differences  are  tabulated  in  terms  of:  root  mean  square  coefficient 

differences,  percentage  differences,  root  mean  square  undulation,  and  root 
mean  square  anomaly  differences.  For  the  degree  8  solution  the  maximum 
undulation  difference  is  67cm  at  degree  8  and  for  the  12  degree  solution  it  is 
91cm  at  degree  11.  The  coefficients  between  9  and  12  of  the  degree  8  solution 
disagree  with  the  GEM9  coefficients  more  that  those  of  the  12  degree  solution. 
Beyond  that  the  differences  are  about  the  same  for  both  solutions. 

The  degree  8  and  12  solutions  in  terms  of  the  implied  coefficients  show 
maximum  differences  at  degree  9  through  14;  beyond  that  the  differences  are 
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small.  The  percentage  differences  are:  at  degree  8,  3%;  at  degree  10,  16%;  at 
degree  12,  26%;  at  degree  14,  11%;  at  degree  16,  8%;  and  at  degree  30,  5%. 
This  shows  that  the  higher  degree  terms  are  not  sensitive  to  the  highest 
degree  used  in  the  adjustment. 

The  adjusted  anomalies  of  the  degree  12  solution  were  developed  into 
harmonic  coefficients  to  degree  60.  The  degree  variances  and  the  root  mean 
square  coefficient  variations  implied  by  these  coefficients  were  compared  to 
GEM9  adjusted  1*  *  1*  mean  anomalies  (Rapp  1977),  and  to  the  geoid  spectrum 
from  altimetry  by  Wagner  (1978)  Kaula's  rule  (10~5/ia)  gave  variations  too 
large  with  respect  to  the  Rapp  solution.  An  excellent  agreement  was  shown 
with  Wagner's  results  from  the  analysis  of  altimeter  data.  The  author  notes 
that  an  adjustment  to  degree  20  would  improve  the  solution,  especially  for 
satellite  orbit  computation, s  but  it  would  take  too  much  space  and  time  on  the 
OSU  computer.  The  effects  of  the  correction  terms  due  to  the  spherical 
approximation,  neglect  of  the  topography  in  computing  the  mean  anomalies,  and 
the  atmosphere  were  not  considered.  Rapp  considers  as  most  critical  the 
terrain  correction,  which  reached  6%  at  degree  39  in  his  previous  Btudy  of  5* 
anomalies  solution  (Rapp  1977). 

The  12  degree  adjustments'  mean  anomalies  were  converted  later  to 
potential  harmonic  coefficients  to  degree  180  using  an  efficient  algorithm  of 
Rizos  (1979)  as  discussed  by  Rapp  (1979a).  The  deficiency  of  this  180  *  180 
set  of  coefficients  was  that  the  coefficients  above  12  degrees  (from  GEM9) 
were  omitted  from  the  adjustment,  therefore  making  it  inadequate  for  orbit 
computations.  Studies  and  experimental  computations  were  carried  out  (Rapp 
1980)  to  facilitate  the  incorporation  of  higher  degree  coefficients  without 
significant  additional  computer  requirements.  An  approximate  procedure  using 
certain  assumptions  was  used  to  generate  a  set  of  180  *  180  coefficients;  this 
was  compared  to  the  rigorous  solution  (Rapp  1978).  The  average  percentage 
difference  was  8.6%,  the  root  mean  Bquare  undulation  difference  was  *  80cm 


and  the  root  mean  square  anomaly  difference  was  2  mgal.  The  percentage 
difference  was  6.0%  for  degrees  2  through  12,  15.4%  for  25  through  36,  and  8% 
for  higher  degrees.  Relative  to  the  data  noise  percentage  error  (about  50%) 
the  difference  between  the  rigorous  and  the  approximate  solutions  was 
considered  very  small.  The  approximate  combination  method  was  used  for 
Rapp's  following  model  (Rapp  1981). 
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This  model  used  the  most  recent  available  input  material  in  an  approximate 
type  of  adjustment.  A  selected  set  of  harmonic  coefficients  up  to  order  and 
degree  36  were  combined  with  the  new  64800  1*  *  1*  surface  mean  anomalies. 
The  results  were  the  adjusted  1*  *  1*  anomalies  and  the  adjusted  coefficients 
corresponding  to  the  a  priori  coefficients.  The  adjusted  anomalies  were 
developed  into  spherical  harmonic  coefficients  through  order  and  degree  180. 

The  satellite  observations  are  represented  by  a  set  of  harmonic 
coefficients  merged  from  selected  solutions.  For  these  "a  priori"  coefficients 
all  available  solutions  were  considered.  These  models  are  briefly  described  in 
Rapp's  report  (1981).  Likewise  the  considered  zonal  coefficients  and  resonance 
terms  are  identified  in  the  report.  Before  selection  and  merging  of  the 
coefficients,  1*  *  1*  mean  anomalies  were  computed  from  each  set  and 
compared  with  a  terrestrial  and  GEOS-3  1*  »  1*  data  set.  The  differences 
between  the  terrestrial  and  harmonic  coefficient  values  are  due  mainly  to  the 
difference  in  the  frequency  constraint  of  the  two  types  of  data.  The  small 
differences  between  the  various  harmonic  solutions  was  considered  not  to  be 
significant.  The  merge  procedure  was  to  form  a  weighted  average  for  a 
coefficient  from  the  PGS  S-2,  PGS  1331,  PGSL-1,  and  the  "miscellaneous” 
coefficient.  The  weighting  of  the  PGSL-1  and  miscellaneous  coefficients  was 
done  by  the  use  of  the  standard  deviations  for  each  coefficient.  For  the 
other  two  sets  the  standard  deviations  were  taken  to  be  0.9  of  the 
corresponding  GEM9  coefficient.  If  it  did  not  exist,  the  standard  deviation  of 
the  miscellaneous  coefficient  was  used. 

The  final  set  of  coefficients,  called  "SET  1",  was  checked  by  comparing 
the  computed  1*  *  1*  anomalies  with  the  corresponding  surface  data.  Rapp’s 
(1981)  report  shows  the  mean  square  anomaly  differences  between  GEM9,  PGS 
S-2,  PGS  S-4,  PGS  1331,  PGSL  1,  and  "SET  1"  harmonic  coefficient  fields  and 
the  Terrestrial/GEOS-3  Data  Set.  There  is  no  significant  change.  There  are 
also  no  significant  differences  between  the  anomaly  degree  variances  for  the 
GEM9  and  "SET  1"  coefficients  (Table  3  of  Rapp's  1981  report). 

The  1*  «  1*  mean  gravity  anomalies  used  in  the  combination  were  obtained 
by  merging  an  updated  set  of  land  1*  *  1*  anomalies  (42585  values) 
and  1*  ■  1*  anomalies  computed  from  adjusted  SEASAT  altimeter  data 
(Rowlands  1981).  The  merged  set  contained  56751  anomalies.  All  anomalies 
were  referred  to  the  1980  Geodetic  Reference  System  (Moritz  1980).  In 


comparing  the  new  set  to  the  previous  set  merged  with  G  EOS-3  mean 
anomalies,  a  root  mean  square  difference  of  *7.5  mgal  was  found  on  the  basis 
of  52972  common  values.  The  standard  deviation  of  all  values  in  the  new  set 
is  *10  mgal  (Rapp  1981). 

In  the  application  of  equations  relating  coefficients  and  anomalies  or  the 
quadrature  formula  used  for  the  estimation  of  the  potential  coefficients,  64800 
1*  *  1*  anomalies  are  required.  For  the  missing  8049,  values  the  mean 
anomalies  implied  by  the  coefficients  to  degree  36  of  the  model  "SET  1"  were 
used.  Considering  the  standard  deviation  of  the  implied  anomalies  to  be  *30 
mgal,  the  root  mean  square  standard  deviation  of  the  complete  aet  is  about 
*15  mgal.  Rapp  used  *20  mgal,  allowing  for  the  errors  in  the  anomalies 
computed  by  geophysical  correlations  in  the  unsurveyed  areas  (the  adjustment 
technique  used  requires  that  the  accuracy  figure  for  all  blocks  be  assumed 
the  same). 

The  process  of  combination:  If  a  set  of  mean  anomalies  is  given,  the 
coefficients  and  anomalies  can  be  related  by  a  quadrature  formula  derived  by 
Colombo  (1980).  He  also  showed  that  the  formula  gives  results  almost  as  good 
as  those  of  more  complicated  optimal  estimation  techniques  for  the  harmonic 
coefficients  (Rapp  1981).  The  combination  procedure  used  by  Rapp  was,  in 
essence,  the  computation  of  a  weighted  average  of  the  starting  coefficients 
and  those  implied  by  the  anomaly  data.  The  advantage  of  this  combination 
technique  is  that  normal  equations  are  not  required;  however,  simplifications 
were  made  relative  to  the  observation  weights  for  the  anomalies  in  order  to 
make  a  solution  to  degree  36.  It  was  necessary  to  assume  that  all  the  1*  *  1* 
mean  anomalies  have  the  same  accuracy.  The  principle  of  the  adjustment  was 
initiated  by  Kaula  (1966),  and  the  adjustment  equations  are  given  in  Rapp's 
reports  (1978  and  1981). 

Results:  The  first  results  of  the  combination  were  the  1*  *  1*  mean 

anomalies  and  the  adjusted  coefficients  up  to  order  and  degree  36.  The 
results  were  compared  with  other  recent  models  and  evaluated  to  assess  the 
quality  and  the  improvements  of  the  new  model.  A  summary  of  these 
assessments  are  given  as  follows:  (1)  The  accuracy  of  the  geoid  undulation 
by  degree,  implied  by  the  accuracy  of  the  adjusted  coefficients  from  degree  2 
through  36,  was  *  152cm  for  the  starting  coefficients  and  *87cm  for  the 
adjusted  set,  (2)  the  differences  between  the  a  priori  and  the  adjusted 
coefficients  are  tabulated  in  the  form  of  average  percentage,  root  mean  square 


undulation  and  anomaly  differences  by  degrees.  The  percentage  differences 
are  largest  at  the  higher  degrees  (60%  to  80%  between  degrees  29  to  36),  (3) 
the  anomaly  degree  variances  are  somewhat  higher  at  the  higher  degrees  of 
the  new  solutions. 

The  adjusted  1*  *  1*  anomalies  were  compared  to  the  anomalies  of  Rapp's 
1978  solution.  The  new  -  old  mean  anomaly  difference  was  0.5  mgal,  the  RMS 
difference  was  *11  mgal,  and  the  maximum  difference  was  215  mgal.  The 
absolute  value  of  the  residuals  of  the  adjusted  anomalies  ranges  between  0 
and  18  mgal;  30884  of  the  64800  1*  *  1*  blocks  are  between  0  and  2  mgal  with 
only  3827  above  7  mgal.  The  large  residuals  are  generally  correlated  with  the 
locations  where  the  1*  *  1*  mean  values  have  been  geophysically  predicted. 

The  harmonic  coefficients  were  compared  to  those  of:  Rapp  (1978),  GEM9, 
GEM10B,  GEM10C  and  GRIM3.  The  percentage  difference  between  the  Rapp  1978 
and  1981  models  is  10%  at  degree  7  and  23%  at  degree  10,  increasing  gradually 
to  60%  at  degree  180.  The  difference  between  Rapp's  1981  and  GEM10C  (the 
other  expansion  to  degree  180  published  by  Lerch  and  al.,  in  1981),  is  7%  at 
degree  7,  14%  at  degree  10,  and  120%  near  degree  180.  A  large  difference  in 
the  RMS  anomaly  of  *7.3  mgal  is  tabulated  between  these  two  180*  *  180* 
solutions.  The  anomaly  differences  between  Rapp's  (1981)  model  and  the  other 
compared  expansions  range  from  *3.6  to  *9.1  mgal. 

The  standard  deviation  of  the  adjusted  coefficients  (order  and  degree  36) 
is  obtained  in  the  adjustment;  for  the  higher  degree  coefficients  (not  part  of 
the  adjustment),  accuracy  estimates  were  assigned  composed  from  data  noise 
(*20  mgal)  and  sampling  error  due  to  the  finite  size  of  the  anomaly  blocks.  At 
degree  50  the  percent  error  due  to  data  noise  is  52%  while  the  sampling  error 
is  only  4%;  at  degree  175  the  data  noise  part  of  the  percent  error  is  138%  and 
the  sampling  error  part  is  52%. 

Geoid  undulation  accuracy  and  gravity  anomaly  accuracy  by  degree  and 
cumulativity  are  plotted  on  graphs  in  the  report.  These  are  compared  to 
GEM9  and  to  the  a  priori  coefficients  respectively.  The  percentage  error  of 
the  adjusted  coefficients  are  shown  on  a  graph  with  the  same  errors  of  GEM9 
and  the  a  priori  coefficients.  This  graph  Bhows  that  the  Rapp  1981  solution 
reaches  100%  error  at  degree  120;  the  a  priori  coefficients  do  so  at  about 
degree  30  and  GEM9  at  about  degree  20. 

Anomaly  degree  variances  are  plotted  for  GEM10C,  GRIM3,  and  Rapp  81 
adjusted  coefficients  and  those  implied  by  Kaula's  rule.  GRIM3  values  are 


higher  than  the  other  two  after  degree  20.  The  summations  of  the  degree 
variances  from  degree  2  to  36  for  each  model  are  the  following:  GEM  10C,  239 
mgal;  GRIM3,  265  mgal;  Rapp  81,  228  mgal.  After  about  degree  50,  the  degree 
variances  of  the  GEM  10C  expansion  are  lower  than  Rapp’s  solution  by  about  1 
mgal. 

For  testing  the  model  for  orbital  computations,  a  number  of  orbital 
comparisons  were  made  at  Goddard  Space  Flight  Center  by  Frank  Lerch  and 
James  Marsh.  The  results  of  these  tests  indicate  that  Rapp’s  model  performed 
comparably  with  some  solutions  in  existence  at  that  time,  such  as  GEM9, 
GEM10B,  and  PGSL  1.  However  the  author  points  out  that  additional  testing  is 
needed  to  obtain  a  better  picture  of  the  performance  of  this  model  in  orbital 
work.  It  is  also  stated  in  the  report  that  the  intent  was  to  produce  a  general 
field  and  not  one  specifically  tailored  to  a  satellite. 

In  the  conclusions  are  mentioned  the  required  improvements  generally 
given  also  by  other  authors  such  as:  surface  topography  corrections  for  the 
altimetry  data;  better  mean  anomalies  for  the  replacement  of  geophysically 
predicted  values;  and  more  rigorous  combination  procedures. 

2.3  The  Smithsonian  As  trophy  si  cal  Observatory  (SAO  Gaposchkin)  Models 

2.3.1  SAO-79  model 

This  is  a  description  of  the  earth’s  gravity  field  in  spherical  harmonics 
through  degree  and  order  30  (Gaposchkin  1980).  Three  types  of  data  have 
been  used  in  the  combination:  satellite  tracking,  terrestrial  gravity 

observations,  and  satellite-altimeter  data. 

Satellite  laser-ranging  data  on  10  satellites  from  15  globally  distributed 
stations  were  used.  The  accuracy  of  the  observations  ranges  from  5m  (1971) 
to  5cm  for  the  1975  data.  The  normal  equations  from  a  previous  model 
(Gaposchkin  and  Mendes  1977)  involving  the  same  satellites  and  observations 
were  used  in  the  adjustment  of  this  solution.  The  procedure  for  orbit 
analysis  is  given  in  Gaposchkin  (1979a) 

The  1*  *  1*  mean  gravity  anomaly  data  was  obtained  from  the  Defense 
Mapping  Agency,  Aerospace  Center.  These  data  were  merged  with  existing  SAO 
data.  In  total,  1504  550km  -  550km  block  anomalies  were  computed  by 

collocation. 


The  base  for  altimeter  data  was  2116  tracks  with  442,411  data  points.  The 

processed  data  was  sorted  into  1*  *  1*  averages  by  linear  regression  and  the 

RMS  of  the  residuals  were  computed.  The  bad  data  with  a  residual  greater 
than  3a  were  removed.  This  process  removed  65  tracks  and  15,762  data 
points.  The  average  a  of  the  remaining  25,295  1*  *  1*  geoid  heights  was 
2.1m.  The  1*  *  1*  geoid  heights  were  then  formed  into  550km  *  550km  block 
values  by  collocation,  resulting  in  1218  values. 

About  seven  solutions  with  different  weights  were  made,  and  the  one  that 

gave  the  best  agreement  with  all  data  types  was  adopted.  Increasing  the 

influence  of  one  particular  type  of  data  usually  weakens  the  agreement  with 
the  others;  therefore,  the  determination  of  the  relative  weights  is  a  difficult 
task. 

Detailed  comparisons  of  the  various  solutions  with  surface  gravity  data, 
geoid  heights,  and  satellite  orbits  are  tabulated  in  Gaposchkin’s  report  (1980). 
In  addition  to  Gaposchkin’s  seven  solutions,  the  GEM10  model  by  Lerch  and  al. 
(1977)  is  also  included  in  these  comparisons.  The  harmonic  coefficients  of  the 
selected  solution  are  listed  in  the  report.  It  is  also  concluded  that  the  model 
agrees  satisfactorily  with  the  surface  gravity,  altimetry,  and  satellite 
ephemerides  available  at  the  time  of  the  development.  The  author  also  lists 
deficiencies  of  the  basic  information  used,  such  as:  a)  the  satellite  normal 
equations  date  back  to  1976,  with  "poorer  quality  in  both  accuracy  and 
distribution"  with  respect  to  "significant  new  satellite  laser-ranging  data";  b) 
it  is  stated  that  better  reference  models,  orbit  computation  programs,  and 
station  coordinates  would  improve  the  model;  c)  only  a  fraction  of  the 
altimeter  data  has  been  included,  and  some  time  variations  could  be  reduced 
by  averaging  the  data;  d)  many  of  the  ocean  data  (observed  or  predicted) 
should  be  replaced  by  altimeter  data  and  geophysically  interpreted  data  on 
land  areas  must  be  handled  with  caution. 

2.4  The  European  (GRIM)  Models 

2.4.1  The  GRIM  3  Model  (Ch.  Reigber  and  al,  1983a) 

For  the  development  of  this  model  the  following  data  were  used: 

(1)  Orbital  perturbations  for  22  satellites,  consisting  of:  optical  tracking  data 

on  18  satellites,  laser  ranging  on  9  satellites,  and  Doppler  measurements 


on  one  satellite.  These  measurements  were  combined  in  106  arcs  of  5  to 
25  days,  corresponding  to  1028  days  of  continuous  tracking. 

(2)  Condition  equations  for  zonal  and  resonant  harmonics  of  order  11  to  15, 
derived  from  other  studies. 

(3)  Surface  gravity  data  in  the  form  of  1*  *  1*  mean  free  air  gravity 
anomalies  consisting  of  (a)  25001  values  computed  from  measured  or 
predicted  point  values  and  taken  from  Rapp's  Oct.  1979  file;  (b)  27916 
values  over  ocean  areas  derived  from  GEOS-3  altimetry,  (Rapp's  1980  file). 

The  solution  was  obtained  by  combining  the  normal  equation  systems  and 
introducing  constraints  for  some  station  positions  of  "special  importance" 
(larger  weights).  Absolute  constraint  was  introduced  to  impose  well  known 
short  inter-station  vectors  and  the  connections  of  the  station  network  of 
Doppler  data  (MEDOC)  to  the  laser  and  optical  sites. 

The  results  are  spherical  harmonic  coefficients  up  to  order  and  degree  36 
and  coordinates  of  95  tracking  stations.  A  detailed  gravimetric  geoid  was 
computed  from  the  1*  *  1*  surface  gravity  data  and  the  GRIM3  model  by 
8uperimposition  of  the  geoid  computed  by  the  integration  of  Stocks  formula  to 
the  long  wavelength  geoid  from  the  GRIM3  coefficients. 

The  model  has  been  extensively  compared  to  and  evaluated  with  respect  to 
other  recent  solutions,  such  as  GEM10B,  Rapp  81,  and  SA079.  The  comparisons 
show  that  the  model  is  reasonably  close  to  these  models,  especially  to  SA079, 
in  terms  of  zonal  and  resonant  harmonics  of  order  15.  RMS  coefficient 
differences,  anomaly  differences,  and  undulation  differences  by  degree  are 
given  in  plots  between  GRIM3  and  each  of  GEM10B,  SA079,  and  Rapp  81 
solutions.  The  comparisons  show  larger  differences  in  the  long  wavelength 
geoid  (up  to  1  :  10);  the  anomaly  differences  by  degree  increase  with  J.  The 
authors  state  that:  "  there  is  no  evidence  that  one  model  is  better  than  the 
other". 

The  long  wavelength  differences  were  analyzed  by  comparing  potential 
differences  at  high  altitude.  This  analysis  concludes  that  terms  of  degree  2 
and  3  of  GRIM3  are  mostly  involved.  The  authors  suspect  that  the  low  weight 
given  to  the  predicted  gravity  data  is  one  of  the  reasons,  in  addition  to  a 
lack  of  enough  high  and  well  observed  satellite  data  (such  as  LAGEOS).  The 
geoid  heights  were  compared  with  those  of  GEM10B.  The  results  show  the 
following  differences:  maximum,  14.7m;  minimum,  15.8m;  mean,  0.03m;  RMS,  3.7m. 
The  extreme  differences  occur  where  only  predicted  data  was  available.  The 


16 


differences  over  ocean  areas  are  small  (RMS  =  1.5m).  Both  solutions  depend  in 
these  areas  on  GEOS  3  data,  but  they  also  indicate  the  good  quality  of  the 
anomalies  converted  from  GEOS-3  geoid  heights. 

The  geoid  height  over  ocean  areas  were  compared  with  SEASAT 
independent  geoid  heights  along  selected  profiles.  For  these  comparisons, 
profiles  from  GEM10B  and  SA079  were  also  used.  The  conclusion  of  thiB  test 

is  that  GEM10B  fits  the  SEASAT  geoid  better  than  GRIM3  and  that  the 

differences  have  long  wavelength  components  due  to  the  error  in  the  low 

degree  harmonics  previously  discussed. 

For  orbit  computation  GRIM3  solution  with  GEM10B  model  using  GRIM3  and 
Tapley  (1980)  station  coordinates  respectively.  GEOS  1,  GEOS  3,  LAGEOS, 

BEACON  3,  and  Starlette  were  used.  The  conclusions  of  these  tests  are  that 
orbital  errors  are  quite  large  with  both  models;  radial  and  cross-track  errors 
are  very  close,  but  errors  in  the  along-track  component  are  smaller  with  the 
GRIM3  model. 

The  authors  conclude  that  the  model  is  satisfactory  globally,  although  they 
recognize  deficiencies  -  due  to  poor  or  missing  gravity  data  -  of  some  low 
degree  and  order  terms. 


2-4.2  GRIM  3B  Model  (Ch.  Reigber  and  al.  1983b) 

The  data  used  for  the  original  solution  of  GRIM3  was  supplemented  with  the 
following  additional  data: 

(a)  16  months  of  LAGEOS  laser  tracking  data 

(b)  1*  x  1*  SEASAT  altimeter  data  (Rapp) 

(c)  1*  x  1*  land  gravity  anomalies  from  Rapp's  1983  file 

(d)  recent  zonal  and  resonant  harmonics. 

The  new  solution  resulted  in  a  complete  set  of  spherical  harmonic  coefficients 
to  degree  and  order  36,  and  the  coordinates  of  109  tracking  sites. 

A  total  of  151  cures  of  laser  and  optical  measurements  to  21  satellites  were 
used  in  this  solution.  63  arcs  were  LAGEOS  arcs  observed  from  20  laser 
tracking  stations  between  January  1980  and  June  1981.  For  the  LAGEOS  orbit 
analysis  in  addition  to  the  earth’s  gravity  and  lunar-solar  attractions  models 
for  earth  tides,  ocean  tides,  the  earth  albedo  and  along  track  accelerations 
were  used.  For  the  coefficients  all  partials  were  computed  according  to  the 
sensitivity  of  a  particular  orbit  to  perturbations  related  to  specific 
coefficients.  For  LAGEOS  all  partials  for  harmonics  of  degree  and  order  16 


were  computed,  and  for  the  zonals  up  to  degree  19. 

Observation  equations  for  even  and  odd  zonal  harmonics  as  derived  by 
various  authors  were  included  in  the  solution.  In  total,  112  condition 
equations  were  used  for  zonal  harmonics.  For  resonant  harmonics  of  order  11 
to  15,  236  observation  equations  were  added  as  derived  by  various  authors 
between  1974-1981. 

As  was  mentioned  previously,  the  surface  gravity  data  was  contributed  by 
Rapp,  and  consisted  of  18504  land  mean  anomalies  for  1*  *  1*  blocks  and 
38,345  1*  x  1*  sea  anomalies  determined  from  SEASAT  altimetry  covering 
oceanic  and  coastal  areas  between  72*N  and  72*S  latitudes.  The  general 
computation  method  used  for  the  GRIM  models  is  described  in  detail  in  Reigber 
and  al.  (1983a).  Essentially,  it  is  the  combination  of  separate  normal  equation 
systems  derived  from:  a)  satellite  orbit  perturbations  computed  from  tracking 
data  of  close  earth  satellites;  (b)  observation  equations  for  zonal  and  resonant 
terms  obtained  from  analysis  of  secular  and  long  period  perturbations  of 
various  satellites;  and  c)  surface  gravity  data  in  form  of  1*  *  1*  mean 

anomalies  of  both  land  measurements  and  ocean  surface  altimetry  observations. 
The  normal  equations  were  computed  by  directly  using  the  1*  x  1*  values  and 
their  respective  uncertainties  (5  -  15  mgal). 

The  results  of  the  solution  are:  a)  fully  normalized  harmonic  coefficients 
of  the  potential  complete  to  degree  and  order  36;  b)  the  geoid  computed  form 
the  obtained  coefficients;  and  c)  the  coordinates  of  109  tracking  stations. 

The  harmonic  coefficients  were  compared  with  the  coefficients  of  the 
GEM10B,  GEM  L-2,  Rapp  81,  and  SA079  models;  the  RMS  coefficient,  undulation, 
and  anomaly  differences  by  degree  are  plotted  on  graphs.  Comparisons  were 
made  also  with  observed  land  gravity  anomalies,  SEASAT  altimetric  geoid, 
gravity  derived  from  satellite  to  satellite  tracking,  and  for  satellite  orbit 
determinations. 

From  the  detailed  evaluations  and  comparisons  given  in  the  paper,  the 
authors  conclude  that  due  to  the  large  quantity  of  precise  LAGEOS 
observations,  the  longest  wavelength  part  of  the  potential  (harmonic 
coefficients  up  to  order  and  degree  7)  has  been  considerably  improved.  A 
better  fit  to  observed  longitude  accelerations  of  24  hour  satellites  and  to  the 
SEASAT  geoid  has  been  obtained  as  compared  to  GRIM3.  The  utilization  of 
GRIM3B  improved  LAGEOS  orbits;  however,  for  Starlette  and  GEOS  3  the  orbits 
are  better  with  GRIM3  than  GRIM3B.  A  small  improvement  iB  seen  in  the  fit  to 
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surface  gravity  data.  The  coordinates  of  the  109  tracking  sites  range  in 
accuracy  from  10  meters  to  several  centimeters.  The  large  errors  are  for  the 
old  optical  stations  and  the  most  accurate  are  the  coordinates  of  the  20 
LAGEOS  laser  tracking  sites.  These  are  compared  with  the  GEM  L-2  set.  No 
major  difference  was  found  between  the  two  sets  of  coordinates. 

The  authors  feel  that  the  longest  wavelength  part  of  the  model  is  "reliably 
determined".  As  far  as  the  rest  of  the  spectrum  is  concerned,  improvements 
are  necessary. 

2.4.3  The  GRIM3-L1  and  the  GRIM3-MIP  models 

Since  the  derivation  of  the  GRIM3B  model,  two  additional  models  were 
derived  by  the  same  group:  GRIM3-L1  (Reigber  and  al.  1984),  and  GRIM3-MIP 
(Reigber  and  al.  1984a).  These  models  are  essentially  recombinations  of  the 
observation  material  used  for  GRIM3B,  with  different  weighting  of  particular 
data  sets  according  to  the  dedication  of  the  model.  GRIM3-MIP  (preliminary 
model)  contains  Doppler  tracking  data  of  TRASIT  (three  months  data  from  the 
MEDOC-1  tracking  campaign).  This  model  was  tailored  for  the  processing  of 
MEDOC-2  Doppler  data.  Although  the  description  of  GRIM3-L1  is  still  not 
available  in  the  open  literature,  in  Reigber  and  al.(2984)  it  is  stated  that  this 
model  has  a  well-balanced  weighting  and  it  is  not  tailored  to  a  particular 
satellite.  GRIM3B  (Reigber  and  al.  1983b)  is  considered  a  LAGEOS  model.  The 
latest  three  GRIM  models:  GRIM3B,  GRIM3-L1,  and  GRIM3-MIP,  have  been 
compared  with  GEM10B  and  GEM  L-2  models,  with  observed  gravity  data  and 
for  orbit  determinations.  The  conclusion  is  that  GRIM3-L1  and  GRIM3-MIP  give 
a  better  fit  to  any  data  source  than  previous  solutions  (GRIM  models).  A 
significant  improvement  was  achieved  in  the  determination  of  low  order 
harmonics.  Further  improvements  are  expected  from  more  accurate  and  better 
data  coverage  (Reigber  and  al.  1984a) 

2.5  Special  Models 

The  Geopotential  models  reviewed  in  the  preceding  sections  belong  in  the 
group  of  "general  purpose"  models  intended  for  both  terrestrial  use  and 
satellite  orbit  computations.  Most  of  these  models  were  derived  by  combination 
of  satellite  tracking  and  various  surface  gravity  data.  The  weights  of  the 


contributing  information  were  usually  determined  by  experimental  solutions 
with  the  goal  that  the  resulting  coefficients  satisfy  the  best  all  types  of  input 
data  and  fit  several  satellite  orbits.  An  exception  is  the  GEM9  model,  which 
contains  only  satellite  observations;  however,  this  model  provided  the  satellite 
tracking  data  for  many  subsequent  general  and  special  combination  solutions. 
Therefore,  the  GEM9  and  its  refined  model,  the  GEM  L-2,  are  discussed  in  the 
preceding  group  of  the  general  models.  Special  models  are  developed  for  the 
computation  of  satellite  orbits  or  more  specifically  "tailored"  for  the  orbit  of 
specific  satellites  (Gravity  models  used  for  the  computation  of  the  Navy's 
navigation  satellites,  or  those  developed  for  G EOS-3  and  SEASAT). 

2.5.1  Special  Models  of  the  Goddard  Space  Flight  Center:  PGS-1.  PGS-S2. 

PGS-S3  and  PGS-S4  (Lerch  and  al.  1982b) 

The  purpose  of  this  development  was  to  improve  the  geopotential  model  for 
the  computation  of  SEASAT  ephemerideB  to  achieve  a  10cm  accuracy  for  the 
radial  component  of  the  satellite  position.  The  altimeter  data  precision  of 
SEASAT  is  somewhat  better  than  10cm;  therefore,  the  full  utilization  of  the 
altimetric  precision  for  geodetic  and  oceanographic  work  required  the  10cm 
accuracy  for  the  radial  position  of  the  spacecraft. 

The  ephemerides  computed  with  the  GEM9  model  and  the  tracking  station 
coordinates  by  Marsh  and  al.  (1977)  indicated  an  error  of  3-5m  RMS  in  the 
radial  position.  First,  the  normal  equations  of  GEM9  were  combined  with  eight 
3-day  laser  orbits  of  SEASAT.  This  solution,  called  PGS-S1,  was  complete  to 
degree  and  order  30  with  additional  terms  to  degree  36.  This  solution 
provided  some  improvement,  but  the  radial  errors  remained  several  meters  in 
magnitude.  In  the  second  step  of  development,  the  Unifid  S.  Brand  (USB) 
data  for  the  same  eight  arcs  and  for  one  additional  arc  of  SEASAT  data  was 
added  to  form  PGS-S2  model.  The  radial  accuracy  of  this  model  was  improved; 
however,  large  differences  showed  in  accuracy  between  the  regions  with 
adequate  tracking  station  and  those  with  sparse  tracking  sites  or  none  at  all. 
In  the  Northwest  Atlantic  region  the  accuracy  was  about  lm,  and  in  the  South 
Atlantic  and  Pacific  areas  about  3m.  The  RMS  value  was  1.8m. 

Analyses  indicated  that  tracking  data  accuracy  is  adequate  to  obtain  a 
50cm  radial  accuracy,  provided  a  very  accurate  force  model  is  given  (Lerch 
and  al.  1982).  There  was  no  more  ground  tracking  data  for  the  improvement 
of  the  geopotential  model,  thus  it  was  decided  to  use  GEOS-3  altimeter  data. 


The  GEM10B  model  consisting  of  GEM9  +  GEOS-3  +  5*  *  5*  surface  mean 

anomalies  was  used.  The  normal  equations  of  this  model  were  combined  with 
SEASAT  tracking  normal  equations  into  a  new  model,  PGS-S3,  developed  to 
degree  and  order  36.  The  authors  state  that  "a  quasi-stationary  sea  surface 
topography  model  was  not  used,  therefore  there  may  be  some  aliasing  within  a 
few  low-degree  and  low-order  coefficients  of  the  geopotential  due  to  these 
unmodeled  effects'*.  They  consider  this  effect  significant  for  SEASAT  orbit 
computation  at  the  10cm  level.  It  is  recognized  that  the  global  validity  of  the 
coefficients  may  be  sacrificed,  but  it  is  pointed  out  that  the  purpose  of  the 
development  was  not  to  create  a  general  model,  but  a  specific  one  for  the 
SEASAT  orbit  requirements. 

The  errors  of  the  SEASAT  radial  positions  improved  to  about  1.2m  RMS,  a 
substantial  improvement  over  PGS-S2.  This  model  has  been  used  by  the  Jet 
Propulsion  Laboratory  for  the  computation  of  the  ephemerides  for  the  final  set 
of  the  released  SEASAT  altimeter  data.  The  PGS-S3  coefficients  and  the 
coordinates  for  the  laser  and  USB  stations  are  listed  in  the  report  (Lerch  and 
al.  1982b). 

The  GEOS-3  altimeter  data  has  a  number  of  limitations  as  compared  to  the 
SEASAT  data  (altimeter  precision,  lack  of  data  recording  system  restricting  the 
altimetry  to  shorter  arc  lengths).  The  SEASAT  altimeter,  with  10cm  precision 
and  with  a  data  recording  system,  permitted  continuous  coverage  of  the 
oceans.  A  set  of  9600  globally  distributed  observations  covering  12  days  were 
combined  with  the  data  of  PGS-S3  model.  This  combination  resulted  in  PGS-S4, 
a  set  of  coefficients  (36  x  36)  and  adjusted  station  coordinates.  This  model 
reduced  the  SEASAT  radial  error  to  about  70cm  RMS 

The  basic  technique  for  the  combinations,  like  for  the  other  GEM  models, 
was  a  weighted  least  squares  adjustment.  The  subsequent  solutions  of  this 
set  of  models  have  been  progressively  produced  by  the  addition  of  new  data. 
For  PGS-S1  and  PGS-S2,  where  only  laser  and  USB  data  were  added,  it  was 
necessary  to  use  a  "modified  least  squares  approach".  This  technique  iB 
described  as  "using  the  a  priori  statistics  of  the  gravity  field  to  provide 
stability  for  the  recovery  of  higher  degree  coefficients  and  to  reduce  the 
effect  of  aliasing".  The  method  minimizes  both  the  signal  and  the  noise,  thus 
controlling  the  over-adjustment  which  occurs  in  the  case  of  simple  least 
squares  adjustment,  when  high  correlation  exists  between  high  degree  and 
order  coefficients  in  the  presence  of  noise  (observation  residuals).  The 


mathematics  is  described  in  Lerch  and  al.  (1979),  in  the  discussion  of  GBM9. 
[In  the  paper  "Goddard  Earth  Models  for  Oceanographic  Applications  (GEM10B 
and  IOC,  Lerch  and  al.,  1981)  the  "modified  least  squares"  technique,  called 
“least  squares  collocation",  was  also  used  in  the  adjustment  of  GEM10.] 


The  SEASAT  gravity  models  have  been  tested  extensively  regarding  their 
performance  in  the  determination  of  the  radial  position  of  the  satellite.  One 
group  of  these  tests  was  the  use  of  altimeter  data  at  track  intersections  for 
the  determination  of  radial  position  errors;  the  other  was  the  comparison  with 
the  independent  orbital  computations  for  SEASAT  by  the  Naval  Surface 
Weapons  Center  (NSWC)  from  TRANET/Georeceiver-Doppler  data  (Colquitt  and  al. 
1980).  The  altimeter  crossover  test  gave  1.4m  RMS  for  the  NSWC  "smoothed" 
ephemeris  versus  the  1.1m  RMS  of  PGS-S4. 

2.5.2  The  Rapp  1977  model 

Another  "special  model"  is  the  set  of  potential  coefficients  to  degree  52 
computed  by  Rapp  (1977a)  from  terrestrial  anomalies  only.  The  purpose  of 
this  experiment  was  not  to  produce  the  best  possible  model  but  to  analyze 
certain  aspects  of  the  terrestrial  gravity  field. 

From  38,406  1*  *  1*  mean  free  air  anomalies,  1507  5*  equal  area  anomalies 
were  computed,  with  an  additional  147  predicted  values  -  a  total  of  1654 
globally  distributed  anomalies  were  available.  The  summation  formulae  were 
used  to  derive  the  harmonic  coefficients.  A  smoothing  operator  was  used  in 
these  computations  and  was  found  to  significantly  effect  the  higher  degree 
coefficients.  For  the  determination  of  the  terrain  effect  several  different 
terrain  correction  models  were  UBed.  It  was  indicated  that  the  terrain 
correction  to  low  degree  coefficients  is  on  the  order  of  10X  to  25X.  It  was 
found  that  the  corrected  potential  coefficients  did  not  agree  as  well  as  the 
uncorrected  coefficients  with  the  satellite  derived  GEM7  coefficients.  The 
percentage  accuracy  by  degree:  8.6X  at  degree  3,  53X  at  degree  12,  86X  at 
degree  30,  and  88X  at  degree  52. 

Mean  anomalies  were  computed  from  the  potential  coefficients  and  compared 
to  the  original  anomalies.  The  difference  was  *8.6  mgal  at  degree  16,  *4.7 
mgal  at  degree  36,  and  *2.6  mgal  at  degree  52.  Obviously,  better  coverage 
and  better  quality  gravity  data  would  have  given  better  results. 


ial  Models  of  the  Naval  Surface  Weapons  Center 


These  models  consist  of  harmonic  coefficients  computed  through  nineteenth 
order  and  degree  with  some  additional  terms  or  deletions  in  different 
variations.  The  models  are  optimized  for  the  polar  orbits  of  the  Navy 
Navigation  Satellites  and  generally  were  not  intended  for  use  for  other 
satellites  and/or  terrestrial  purposes.  The  principal  geodetic  ubo  of  these 
models  was  for  the  computation  of  Doppler  satellite  ephemerides,  and  they 
were  revised  many  times  from  the  original  inception  of  the  system.  Between 
1967  and  1972,  the  gravity  fields  designated  as  NWL-8D,  NWL-8H,  and  NWL-9B 
have  been  used  for  the  Doppler  system.  Some  information  regarding  these 
models  can  be  found  in  (Anderle  and  al.  1976).  The  gravity  fields  used  since 
1972  are  listed  below  with  information  available  in  the  open  literature.  Some 
additional  information  was  obtained  from  Richard  J.  Anderle  of  Naval  Surface 
Weapons  Center  (NSWC). 

2.5.3. 1  World  Geodetic  System  1972  (WGS72)  Gravity  Model.  This  model  was 
derived  simultaneously  with  the  development  of  the  World  Geodetic  System  1972 
(WGS72)  by  a  Department  of  Defense  (DOD)  working  group  (Seppelin  1974). 

The  model  consists  of  harmonic  coefficients  complete  through  degree  and 
order  19,  with  additional  zonal  coefficients  through  degree  24  and  resonance 
terms  through  order  27.  The  data  used  for  the  development  was  that  of 
Doppler  observations  of  the  Navy  Navigation  Satellites  (NAVSATs)  taken  at  the 
semi-permanent  Doppler  network  (TRANET  stations).  This  is  about  85%  of  the 
Doppler  data;  15%  was  observed  by  equipment  in  mobile  vans  at  over  120 
worldwide  distributed  locations;  optical  observations,  laser  data  and  surface 
mean  anomalies  in  the  form  of  normal  equations  of  the  Smithsonian  Earth  Model 
II,  derived  by  the  Smithsonian  AstrophyBical  Observatory  (SAO);  and  surface 
gravity  data  in  the  form  of  410  10*  x  10*  equal  area  mean  free  air  gravity 
anomalies.  Only  45X  of  these  anomalies  were  computed  from  observed  data. 
The  remaining  55X  were  derived  by  "gravity-geophysical  correlation 
techniques". 

For  each  data  set,  a  normal  equation  matrix  was  formed,  then  the  normal 
equation  matrices  were  combined  and  the  resultant  matrix  solved. 

There  is  no  evaluation  or  comparison  of  this  model  with  other  models  in 
the  open  literature. 


2. 5. 3. 2  NWL10-E  Model.  This  model  was  introduced  for  computation  of  all 
precise  NAVSAT  ephemerides  in  January  1973,  replacing  the  previous  gravity 
field  designated  at  NWL9-B.  It  is  based  on  satellite  observations  used  for  the 
development  of  WGS72,  and  gravity  coefficients  are  computed  up  to  order  and 
degree  (n,m)  =  (19,19),  but  selected  terms  below  19th  degree  were  suppressed 
and  higher  order  coefficients  to  (n,m)  =  (29,27)  were  added  (Anderle  1984). 
This  model  produces  the  same  NAVSAT  orbit  as  WGS-72  to  an  accuracy  of  lm. 

2. 5. 3. 3  NWL10E-1.  This  model  replaced  the  previous  NWL10-E  model  in  June 
1977.  It  is  the  same  as  NWL10-E  except  for  the  modification  of  two  27th  order 
resonance  coefficients. 

2.5.3.4  NWL-IG6.  This  gravity  field  waB  a  revision  of  two  pairs  of  coefficients 

of  each  order  from  the  NWL10-E  model  to  obtain  orbits  which  fit  GEOS-3 
Doppler  observations.  It  was  expected  that  comparison  of  altimeter 

measurements  on  intersecting  tracks  will  reduce  the  bias  to  70cm  below  the 
precision  of  the  altimeter  (Anderle  and  al.  1976).  This  set  of  coefficients  was 
used  for  a  year  or  two  for  GEOS-3  and  was  subsequently  replaced.  All 
previous  ephemerides  were  recomputed  UBing  a  Goddard  model,  GEM9  (Anderle 
1984). 

NWL10E1:  This  model  is  the  NWL10E  model  with  improved  14th  and  15th 
order  resonance  coefficients.  ThiB  model  was  adopted  for  precise 
computations  of  SEASAT  orbit  by  NSWC  because  of  a  slight  improvement  in 
radial  residuals  over  those  obtained  with  the  GEM10  model  (Anderle  1984). 

2.5.3.5  World  Geodetic  System  1984  (WGS84).  This  gravity  model  is  part  of  the 
new  geodetic  system  derived  by  the  Department  of  Defense  to  replace  the 
WGS72.  This  model  is  a  combination  adjustment  of  Doppler  tracking  data  of 
near-earth  satellites,  laser-ranging  data  on  LAGEOS  satellite  altimetry  data 
over  the  oceans,  and  mean  gravity  anomalies  of  3*  «  3*  derived  from  surface 
observations.  The  harmonic  coefficients  complete  to  degree  and  order  are 
given  in  Macomber  (1985).  The  model,  however,  is  complete  through  degree 
and  order  42  with  some  higher  degree  terms  (Anderle  1984).  A  Technical 
Report  describing  the  WGS84  is  in  preparation  and  will  be  published  in  1985. 
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3.  PROGRAMS  AND  TECHNIQUES  FOR  THE  IMPROVEMENT  OF  THE  GRAVITY  FIELD 

MODELS 

The  analysis  of  satellite  orbital  perturbations  and  satellite  altimetry 
substantially  advanced  the  description  of  the  earth’B  gravity  field  in  terms  of 
spherical  harmonic  series  and  oceanic  geold  undulations  and/or  mean 
anomalies.  These  new  types  of  data  covering  large  areas  combined  with 
surface  measurements  over  land  areas  facilitated  the  derivation  of  many 
different  sets  of  harmonic  coefficients  expanded  to  various  degrees  according 
to  the  data  used  and  the  purpose  of  the  model.  Some  of  these  models  are 
discussed  in  the  previous  sections. 

One  of  the  principal  error  sources  in  satellite  orbit  computations,  and  in 
turn  in  satellite  altimetry,  is  the  gravity  model.  An  illustration  of  the  need 
for  more  accurate  gravity  models  is  the  accuracy  of  the  current  altimetry. 
The  SEASAT  altimeter  instrumental  precision  is  about  10cm,  on  the  other  hand, 
the  best  radial  accuracy  of  the  orbit  computed  with  the  best  fitting  gravity 
model  is  70cm.  This  is  largely  due  to  the  uncertainties  in  the  gravity  model 
(Lerch  and  al.  1982a). 

Years  ago  the  large  gaps  over  the  oceans  (gravimetrically  empty  areas) 
hindered  the  derivation  of  global  gravity  models  with  meaningful  accuracy; 
today  the  gaps  are  -  excluding  polar  regions  -  over  land  areas  such  as:  (1) 
physically  inaccessible  areas  (high  mountains  and  deserts)  in  Africa,  parts  of 
Asia,  Alaska,  and  South  America;  (2)  countries  where  gravity  data  are 
classified,  even  mean  values  as  large  as  1*  *  1*:  USSR,  Eastern  European 
countries,  and  China  are  examples  (Balmino  1983,  Rapp  1985). 

The  latest  and  the  best  published  1*  *  1*  mean  free  anomaly  coverage  is 
the  "Rapp  January  1983  1*  *  1*  Data  Set"  (Rapp  1983).  Prof.  Rapp  and  his 
staff  collected,  analyzed,  and  processed  gravity  data  during  the  past  ten 
years.  The  data  include:  direct  observations  over  land  and  ocean  areas, 
predicted  values  by  topographic-isostatic  techniques  and  other  geophysical 
correlations. 

The  Jan  83  Rapp  file  contains  44,513  1*  *  1*  mean  anomalies  on  land  and 
sea  derived  from  direct  observations  or  predicted.  These  data  merged  with 
approximately  38000  1*  *  1*  anomalies  estimated  from  SEASAT  altimeter  data 
resulted  in  a  combined  set  of  56000-57000  1*  *  1*  values  used  in  recent 
harmonic  expansions. 


There  is  no  hope  that  the  policy  might  change  in  the  areas  where  gravity 
data  is  classified.  The  physically  inaccessible  areas  may  be  surveyed  by 
airborne  gradiometry  and/or  by  aircraft  carrying  accelerometers  of 
appropriate  design  and  GPS  receivers  for  position  and  aircraft  acceleration 
measurements.  In  other  areas  terrestrial  gravity  data  is  likely  to  be  available 
as  local  or  regional  surveys  are  progressing;  however,  it  cannot  be  expected 
that  large  amounts  of  data  substantially  assisting  global  models  will  be 
produced  by  direct  ground/surface  techniques. 

Substantial  improvements  of  the  gravity  models  can  be  expected  from 
additional  gravity  data  delivered  by  satellite  programs.  Several  programs  and 
various  sensor  instruments  are  under  development  with  very  good  potentials. 
These  programs  are  briefly  reviewed  in  the  sequel. 

3.1  Satellite  Altimetry 

The  sea  surface  information  provided  by  satellite  altimetry  contributed 
dramatically  to  the  knowledge  of  the  geoid,  the  gravity  anomalies  over  the 
oceans  and  in  turn  to  the  improvement  of  the  spherical  harmonic  geopotential 
models.  In  addition  to  the  geodetic  aspects,  altimeter  measurements 
substantially  contributed  to  oceanography  in  the  fields  of:  current  detection, 
ocean  tides,  winds  and  waves,  etc.  Through  the  experience  with  the  Skylab, 
GEOS-3,  and  SEASAT,  the  follow  on  hardware  will  permit  about  2cm  precision 
over  a  period  of  several  years  (Marsh  1983). 

Satellite  altimetry  over  the  past  12  to  15  years  has  become  very  rich  in 
literature.  Several  special  issues  have  been  published  devoted  to  the  various 
scientific  fields  related  to  satellite  altimetry.  A  comprehensive  review  and  a 
long  reference  list  is  given  by  J.G.  Marsh  in  the  Reviews  of  Geophysics  and 
Space  Physics  (1983).  From  the  point  of  view  of  geodesy,  the  shape  of  the 
mean  sea  surface  approximating  the  geoid  and  the  recovery  of  mean  gravity 
anomalies  were  the  primary  objectives  of  the  analyses  of  altimeter  data.  An 
estimate  of  the  sea  surface  relative  to  the  geoid  for  the  North  Atlantic  Ocean 
waB  published  by  Winch  (1981).  Similar  work  was  done  by  Stommel  and  al. 
(1978),  and  for  the  global  ocean  surface  by  Levitus  (1982).  Numerous  analyses 
have  been  done  to  obtain  the  sea  surface  heights  (geoid)  and  mean  gravity 
anomalies  from  both  GEOS-3  and  SEASAT  data.  Some  of  these  analyses  are: 


Rapp  (1977,  1979a,  1979b,  1979c,  1982a,  1982b,  1983b),  Cruz  (1983),  Hadgigeorge 
and  al.  (1979,  1981),  Kearsley  (1977),  Kahn  and  al.  (1979a,  1979b),  Marsh,  J.G. 
and  al.  (1980),  Marsh,  J.G.  and  Martin  (1982),  Liang  (1983),  Rowlands  (1981), 
Marsh,  J.G.  and  al.  (1983). 

Satellite  altimeter  data  were  used  in  combination  solutions  for  global 
gravity  models  by:  Rapp  (1978,  1981),  Gaposchkin  (1980),  Hadgigeorge  and 
Blaha  (1981),  Lerch  and  al.  (1981,  1982a),  Reigber  and  al.  (1983a,  1983b,  1984), 
and  others.  These  models  are  reviewed  in  the  preceding  sections  of  this 
report. 

The  most  extensive  gravity  anomaly  work  from  altimetry  data  was 
performed  by  Rapp  and  his  research  associates  at  the  Ohio  State  University. 
Some  results  of  their  work  is  reviewed  below:  The  computation  of  mean 

anomalies  and  geoid  heights  from  GEOS-3  data  is  described  in  Rapp  (1977, 
1979a,  1979b).  The  altimeter  data  were  received  from  NASA  Wallops  Flight 
Center  up  to  Sept.  1977.  After  editing,  the  data  set  contained  419,216  frame 
measurements  in  1976  arcs.  The  orbital  error  and  altimeter  bias  were  removed 
by  an  adjustment  using  a  first  degree  polynomial  in  time  for  the  orbit  and 
bias  error,  and  introducing  crossover  observation  equations  (Rapp  1979b). 
From  a  set  of  adjusted  geoid  undulations  geoid  undulation  maps,  mean  gravity 
anomalies  and  mean  undulations  were  produced  by  least  squares  collocation  in 
12,144  1*  *  1*  blocks  and  for  377  blocks  of  5*.  The  accuracy  of  a  1*  *  1* 
mean  anomaly  was  predicted  as  *7  mgal  and  it  was  *3  mgal  tor  a  5*  block. 
These  figures  represented  an  improvement  factor  of  2  in  relative  accuracies 
obtainable  from  the  GEM9  model. 

Rapp  and  his  research  associates  at  the  Ohio  State  University  (OSU) 
analyzed  SEASAT  altimeter  data  covering  all  ocean  areas  between  70*  North 
and  South  latitudes.  The  analysis  consisted  of:  the  editing  of  the  basic  data 
produced  by  the  Jet  Propulsion  Laboratory  and  provided  to  OSU  by  the 
National  Geodetic  Survey;  the  adjustment  of  the  data  using  the  crossing  arc 
criteria;  and  the  recovery  of  mean  gravity  anomalies  and  sea  surface  heights 
above  the  reference  ellipsoid.  Based  on  a  set  of  criteria  established  from 
previous  experience,  1,085,006  values  were  deleted  from  the  total  of  4,429,491 
values.  A  second  edit  was  performed  by  fitting  the  altimeter  sea  surface 
heights  to  the  geoid  undulations  derived  from  Rapp's  geopotential  model  to 
degree  180  (Rapp  1978).  This  edit  deleted  6420  data  points.  1667  duplicate 
observations  were  found,  so  that  after  editing,  3,336,398  observations  remained 


in  5036  arcs  (Rapp  1982b)  The  adjustment  was  performed  by  Rowlands  (1981) 
using  the  crossing  arc  technique.  Due  to  the  limitations  of  the  computer,  first 
a  primary  set  of  global  arcs  were  adjusted,  then  holding  these  fixed,  four 
regional  areas  followed.  The  average  crossover  difference  before  the 
adjustment  was  *1.5m,  after  the  adjustment  it  was  reduced  to  *28cm.  The 
adjusted  SEASAT  surface  heights,  by  ignoring  the  effects  of  sea  surface 
topography,  can  be  regarded  as  geoid  undulations.  From  the  undulations 
34,973  1*  *  1*  and  1178  5*  *  5*  mean  anomalies  and  geoid  heights  were 
computed  by  least  squares  collocation,  previously  utilized  for  GEOS-3  data  and 
described  in  Rapp  (1979a).  The  comparison  with  common  GEOS-3  values 
resulted  in*7.8  mgal  difference  (standard  deviation  between  27221  common,  1*  * 
1*  anomalies)  for  the  anomalies,  which  corresponds  to  *87cm  for  the 
undulations.  Comparisons  of  the  common  5*  *  5*  values  gave  *2.2  mgal  and 
*76cm  differences  for  the  anomalies  and  undulation,  respectively  (Rapp  1983b). 
Point  sea  surface  heights  computed  from  the  adjusted  data  were  compared 
with  similar  GEOS-3  values.  The  mean  difference  was  1.3m  with  standard 
deviation  of  ‘37cm.  Later  it  was  discovered  that  some  crossover  data  was  not 
used  in  the  adjustment.  The  northeast  Pacific  area  was  primarily  affected  by 
this  error.  It  was  determined  that  the  effect  on  1*  *  1*  mean  anomalies  was 
of  the  order  of  *6  mgal,  and  on  the  surface  heights  *40cm.  Outside  this  area 
the  errors  are  significantly  smaller  (Rapp  1983b).  A  number  of  regional 
adjustments  repeating  the  original  Rowland  (1981)  adjustment  were  carried  out 
in  specific  areas.  Information  on  these  adjustments  are  given  in  Rapp  (1982a) 

Another  work  accomplished  at  OSU  with  GEOS-3  and  SEASAT  altimeter  data 
is  the  "Adjustment  and  Combination  of  GEOS-3  and  SEASAT  Altimeter  Data"  by 
Liang  (1983).  The  previous  GEOS-3  adjustment  of  Rapp  (1979b)  included  only 
3275  arcs.  In  1981  OSU  received  from  NOAA-NGS  the  3.5  year  GEOS-3  data  set 
including  10,520  arcs  from  April  1975  to  December  1978.  This  data  set  iB 
described  in  detail  in  Agreen  (1982).  This  revised  data  set  was  processed  at 
OSU  with  the  goal  that  the  "new"  GEOS-3  data  be  on  the  same  system  with  the 
SEASAT  adjustment  so  that  the  two  sets  of  data  can  be  combined.  A  new 
editing  and  crossover  computing  procedure  was  used  which  has  several 
advantages  versus  the  previously  used  procedure  at  OSU. 

The  combined  set  of  data  is  a  very  dense  coverage  including  all  available 
sea  surface  height  data  from  both  GEOS-3  and  SEASAT  missions.  The 
combined  data  set  also  has  a  very  good  distribution  due  to  the  different 


inclinations  of  the  two  satellites. 


The  differences  between  -he  adjusted  GEOS-3  and  adjusted  SEASAT  1*  * 
1*  surface  heights  in  the  sense  of  GEOS-3  -  SEASAT  is  5cm  with  standard 
deviation  *52cm  (28,810  1*  «  1*  blocks  compared). 

The  comparison  of  the  combined  data  set  with  the  adjusted  GEOS-3  and 
the  adjusted  SEASAT  data  gave  the  following  results:  The  average  1*  *  1* 
mean  surface  height  difference,  GEOS-3  -  combined,  considering  two  different 
areas,  is  3cm  with  an  average  standard  deviation  of  *32cm.  The  average  1*  * 
1*  mean  surface  height  difference,  SEASAT  -  combined,  is  -1cm  with  an 
average  standard  deviation  of  *27cm. 

For  the  areas  of  the  "old"  GEOS-3  coverage  Kearsley  (1977)  prepared  26 
geoid  undulation  maps  from  the  adjusted  data  by  Rapp  (1977).  The  final 
adjusted  data  of  SEASAT  (Rowland  1981,  Rapp  1982a)  were  used  by  Rapp 
(1982b)  to  construct  53  maps  of  sea  surface  heights  at  a  contour  interval  of  2 
meters.  The  heights  are  referred  to  the  1980  Geodetic  Reference  System 
Bllipsoid:  a  =  6378137m,  f  =  1/298.257222.  The  map  is  centered  of  data  on  a 
1*  *  1*  grid  predicted  from  the  adjusted  data  by  a  least  squares  collocation 
procedure. 

The  Air  Force  Geophysical  Laboratory  (AFGL),  formerly  the  Air  Force 
Cambridge  Research  Laboratories  (AFCRL),  developed  a  short  arc  adjustment 
method  for  the  determination  of  the  geoid  and  the  gravity  field  from  satellite 
altimeter  data.  This  method  does  not  require  the  precise  orbit  of  the  altimeter 
satellite  (Brown,  1973,  Blaha  1979,  Hadgigeorge  and  al.  1981).  This  method  was 
used  by  AFCRL  for  the  processing  of  both  GBOS-3  and  SEASAT  altimeter 
observations,  and  it  was  upgraded  to  combine  altimeter  and  gravity  anomaly 
data.  GEOS-3  measurements  were  combined  with  a  set  of  1654  equal  area  5* 
mean  anomalies  (Rapp  1977a),  resulting  in  a  14  *  14  solution  in  spherical 
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harmonics  and  geoid  and  gravity  anomalies  on  a  5*  «  5*  grid  (Hadgigeorge  and 
Blaha  1979).  These  adjustments  (14  *  14)  are  termed  as  "first  phase” 

adjustments  and  they  are  used  as  reference  fields  for  "second  phase" 
solutions  consisting  of  point  maBB  and  collocation  techniques  (Blaha  1984,  Blaha 
and  al.  1984).  The  results  of  the  two  "second  phase”  adjustments  (using  as 
observations  the  residuals  of  the  "first  phase"  14  >  14  adjustment)  represent 
a  2*  resolution,  corresponding  "very  approximately"  to  a  spherical  harmonic 
expansion  of  90  >  90.  The  results  on  a  2*  *  2*  equilateral  grid  were  densified 
into  data  on  a  1*  >  1*  geographical  grid  by  an  "errorless  collocation".  The 
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comparison  of  the  results  for  geoid  undulations  and  gravity  anomalies  was 
accomplished  in  13  large  oceanic  blocks,  containing  12,934  grid  points  (Bessette 
and  Hadgigeorge  1984).  The  r.m.s.  difference  in  geoid  undulation  is  0.45m  and 
the  anomaly  r.m.s.  difference  is  2.8  mgal.  The  average  magnitudes  of  the 
differences  are  33cm  and  2.0  mgal  respectively. 

Efforts  are  in  progress  for  the  improvement  of  the  precision  of  the 
altimeter  system.  These  systems  and  the  expected  improvements  are  described 
in  detail  by  McGoogan  and  al.  (1982),  and  others  in  the  field  of  radar  and 
radio  sciences. 

One  essential  requirement,  to  match  the  altimeter  precision  (currently 
10cm),  is  the  improvement  of  the  radial  orbit  accuracy  of  the  altimeter 
satellite.  This  depends  on  the  improvement  of  current  gravity  potential 
models,  improvement  in  the  global  distribution  and  accuracy  of  tracking  of  the 
satellite,  and  on  the  elimination  of  non-gravitational  force  effects  (air  drag, 
radiation  pressure,  etc.). 

A  number  of  studies  have  performed  and  published  in  a  special  issue  cf 
the  Journal  of  the  Astronautical  Sciences  [28(4)  Oct.-Dec  1980]  on  the  orbital 
accuracy  of  SEASAT.  The  results  can  be  summarized  as  follows:  The  orbital 
errors  of  GEOS-3  have  been  systematically  reduced  for  SEASAT  by  the 
improvements  of  the  geopotential  models,  due  to  satellite  altimetry  data.  The 
several  meters  error  in  GEOS-3  radial  distance  was  reduced  to  75cm  -  1.5m 
r.m.s.  in  the  case  of  SEASAT.  Analyses  by  MarBh,  J.G.  and  Williamson  (1980) 
and  other  studies  also  confirmed  an  apparent  4  meter  difference  in  the  Z 
components  of  tracking  station  coordinates  of  Goddard  Space  Plight  Center  and 
Naval  Surface  Weapons  Center  SEASAT  ephemerides.  The  same  study  indicated 
that  additional  analyses,  in  correlation  with  laser  and  unified  S-band  tracking 
data,  are  expected  to  yield  an  orbital  accuracy  of  50cm  r.m.s. 

There  sure  two  new  altimeter  missions  in  development  in  the  U.S.:  The 
GBOSAT  satellite  of  the  Navy  (Kilgus  and  al.  1981)  and  the  TOPEX  of  NASA. 
The  GEOSAT,  launched  in  1985,  carried  a  SEASAT  type  of  altimeter  with  a 
three  years'  lifetime  and  18  months  of  nominal  mission  time.  The  TOPEX 
mission  according  to  the  plans  will  carry  an  altimeter  with  a  *2cm  precision 
and  will  yield  an  accuracy  of  *14  cm  along  a  measurement  grid  (J.G.  Marsh, 
1973).  A  detailed  analysis  of  the  expected  TOPEX  mission,  contemplated  fro 
1990,  is  given  in  Christensen  (1982).  For  the  tracking  of  TOPEX  and  other 
future  altimeter  missions  the  development  of  a  new  improved  TRANET  Doppler 


network  is  planned.  A  new  tracking  system  is  under  development  at  the  Jet 
Propulsion  Laboratory  (JPL)  using  the  Global  Positioning  system  (GPS).  These, 
together  with  some  improvements  in  the  modeling  of  the  force  fields,  may 

achieve  a  10cm  accuracy  in  the  the  radial  component  of  the  orbits  of  altimeter 

satellites. 

3.2  Satellite- to-Satellite  Tracking  (SST) 

The  technique  of  the  satellite-to-satellite  tracking  for  direct  measurement 
of  the  earth's  gravity  field  originated  in  the  Apollo  program.  Muller  and 

Sjogren  (1968)  used  the  range  rate  tracking  data  of  the  lunar  or  biters  as 

"direct  observations"  of  the  acceleration  of  gravity  along  the  line  of  site 
between  the  tracking  station  on  Barth  and  the  satellites  orbiting  the  moon. 
These  observations  showed  circular  mass  concentrations  ("mascons")  over  the 
flatlanda  of  the  moon  (Colombo  1981b).  The  adaptation  of  this  concept  to  the 
Earth's  gravity  field,  was  the  "high-low"  two  satellite  configuration,  realized 
by  the  "low"  orbiting  GEOS-3  (8400km)  tracked  by  Doppler  from  ATS6 
(40,000km)  geosynchronous  satellite.  Wolf  (1969)  proposed  a  pair  of  satellites 
in  a  "low  -low"  configuration,  where  two  satellites  on  a  low,  circular  polar 
orbit  following  as  close  as  possible,  would  measure  the  range  rate  between  the 
two  spacecraft.  This  arrangement  will  yield  a  high  density  global  set  of 
gravity  data  covering  land  and  sea  without  any  gaps.  The  measurements  will 
provide  the  medium  wavelength  spectrum  of  the  gravity  field,  (higher  degree 
spherical  harmonics)  which  cannot  be  determined  accurately  from  ground 
based  tracking.  The  enormous  potential  of  this  technique  generated  interest 
and  research  to  resolve  theoretical  and  practical  problems.  Studies  were  made 
to  optimize  satellite  arrangements,  develop  data  reduction  methods  and  to 
estimate  achievable  accuracies.  Some  of  these  studies  are:  Schwarz  (1970, 
1972),  Hajela  (1974,  1978,  1979),  Douglas  and  al.  (1980),  Rummel  and  al.  (1976), 
Marsh  and  al.  (1981,  1984).  The  above  Btudies  have  been  concerned  with  the 
recovery  of  medium  and  short  wavelength  gravity  field  in  regional  or  local 
areas,  utilizing  GEOS-3  and  ATS-6  "high  -  low"  observations.  As  the 
advantages  of  the  "low  -  low"  configuration  became  obvious,  more  and  more 
studies  and  error  analyses  were  prepared  on  this  topic. 

Rummel  (1980)  estimated  simultaneously  orbital  parameters  and  the  gravity 


field  in  a  study  of  "low  -  low"  SST  experiment,  assuming  an  "optimal” 
situation:  range  rate  precision  *  10*-6  ms~-l,  satellite  distance  of  250km,  and 
an  orbit  altitude  of  200km.  The  a  posteriori  standard  deviations  were:  *0.9m 
for  point  geoid  heights,  *0.7m  for  geoid  height  differences  and  6  to  7  mgal  for 
1*  *  1*  mean  gravity  anomalies.  These  values  compare  well  with  the  results  of 
GEOS-3  altimetry.  A  study  by  Jekeli  and  Rapp  (1980)  estimated  the  accuracy 
of  mean  anomalies  and  geoid  undulations  for  various  block  sizes  based  on  an 
assumed  mission.  The  accuracy  was  defined  as  a  commission  error  due  to 
measurement  noise  propagation  and  a  truncation  error.  For  a  "low  -  low" 
mission  of  six  months  duration,  at  an  altitude  of  160km  and  range  rate  data 
noise  of  *1  micrometer/sec.  for  four  second  integration  time,  it  was  estimated 
*2.3  mgal  accuracy  for  1*  *  1*  mean  anomalies  and  4.3cm  for  1*  *  1*  mean 
geoid  undulations. 

Colombo  (1981b)  made  an  error  analysis  of  the  global  geopotential  model 
from  a  "low  -  low”  SST  mission.  It  was  assumed  that  two  drag  compensated 
(DISCOS  system)  satellites  were  used.  The  differentiated  range  -  rate  signal 
was  considered  as  equal  to  the  line  of  sight  component  of  the  gravitational 
acceleration.  Both  the  least-squares  adjustment  and  the  least  squares 
collocation  adjustments  were  used.  The  main  results  are  summed  up  by 
Colombo  as  follows: 

IF 

(1)  the  two  satellites  move  on  the  same  polar  circular  orbit  at  160km  altitude, 
at  a  distance  of  300km  from  each  other; 

(2)  the  accuracy  of  the  averaged  range  rate  is  %/~2  *  10~6  ms~l,  the  averaging 
interval  is  4s; 

(3)  residual  data  are  uBed  with  respect  to  a  reference  model  of  specified 
accuracy,  complete  to  degree  and  order  20, 

THEN: 

(1)  The  relative  error  in  the  potential  coefficients  could  be  better  than  IX  up 
to  degree  n  =  130,  better  than  10%  up  to  n  =  210,  and  better  than  50%  up  to 
n  =  270; 

(2)  The  accuracy  of  point  geoid  undulation  implied  by  the  coefficients  could 
be  better  than  5cm  RMS  in  the  band  from  300km  to  40030km  (total  error)  and 
better  than  10cm  in  the  band  from  140km  to  3000km  (total  error). 


It  is  noted,  that  above  degree  200,  the  accuracies  predicted  by  collocation 
are  significantly  better  than  those  by  least  squares.  The  report  considers 
that  the  principles  of  the  study  can  be  applied  to  actual  SST  data,  to  obtain  a 
high  resolution  harmonic  model.  Colombo  recommends  as  adequate  technique 
the  method  in  Colombo  (1981a)  and  that  more  attention  is  to  be  paid  to  global 
data  reduction. 

To  escape  from  the  large  inversion  problems  of  the  least  squares 
collocation,  Kaula  (1983)  developed  an  analytic  scheme  for  inferring  the 
variations  of  the  gravity  field  from  SST.  Each  revolution  is  Fourier  analyzed 
separately,  from  north  pole  to  north  pole  and  east-west  variations  are  inferred 
by  requiring  that  the  potentials  at  the  poles  agree.  The  method  for  data 
analysis  is  outlined  in  the  paper  and  it  is  applied  as  a  test  to  a  pair  of 
satellites  at  160km  and  100km  spacing,  with  5*  data  point  interval  (72  data 
points  per  revolution).  The  assumed  gravity  field  (Gaposchkin  1980)  of 
tesseral  harmonics  to  the  eight  degree  was  completely  recovered  in  three 
iterations  over  64  revolutions.  It  is  demonstrated  that  data  points  at  regular 
intervals  provide  the  opportunity  for  utilizing  techniques  without  massive 
matrix  inversions. 

On  the  basis  of  the  successful  demonstration  of  the  SST  technique  by  the 
results  of  the  GEOS-3/ATS-6  and  of  the  "low  -  low"  simulation  studies  the 
planning  of  gravitational  satellite  missions  started  in  the  mid  1970’s.  NASA  in 
the  US  prepared  the  GRAVSAT  program,  mentioned  earlier  and  the  European 
Space  Agency  (ESA)  came  up  with  a  plan  called  Space  Laser  Low  Orbit  Mission 
(SLALOM)  which  involves  the  simultaneous  tracking  by  laser  interferometry, 
from  the  space  shuttle,  of  two  reflecting  spheres  to  measure  their  relative 
velocities  with  respect  to  the  shuttle  and  to  each  other  (Colombo  1981b) 

The  GRAVSAT  mission  was  reformulated  into  the  NASA  Geopotential 
Research  Mission  (GRM)  with  more  ambitious  accuracy  and  resolution  goals. 
The  GRM  mission  will  be  discussed  later  in  this  report. 

3.3  Satellite  Gradiometry. 

An  approach  with  a  very  good  potential  for  the  determination  of  the  Earth’s 
global  gravity  field  with  higher  resolution  and  accuracy  is  a  gradiometer  of 
high  precision  capable  to  operate  from  a  satellite.  At  least  two  instruments 


with  different  approaches  and  a  number  of  simulation  studies  for  the 
assessment  of  their  performances  are  available.  It  is  expected  that  an 
operational  system  may  be  available  sometime  in  the  1990's. 

One  of  the  instruments  under  development  by  Paik  (1981)  is  the  cryogenic 
gravity  grad  iome  ter.  Superconductivity  and  the  availability  of  SQUID 
(Superconducting  Quantum  Interference  Device)  technology  for  low  noise 
amplifiers  allowed  the  construction  of  a  very  sensitive  and  low  drift  gravity 
gradiometer,  operating  at  liquid  helium  temperatures.  By  employing 
accelerometer  pairs  and  differencing  the  outputs  along  each  coordinate  axis,  a 
"tensor  gravity  gradiometer”  is  obtained,  which  measures  all  the  independent 
components  of  the  gradient  tensor  simultaneously  (five  independent  plus  one 
component).  The  goal  for  the  sensitivity  of  the  instrument  is  10~4E  [IE 
(Eotvos  unit)  =  10“*  S-a] 

Another  satellite-gradiometer  is  under  study  by  a  group  of  French 
scientists,  under  the  name  of  PROJECT  GRADIO.  It  is  described  by  Balmino 
and  al.  (1984).  This  concept  employs  a  gradiometer  composed  of 
micro-accelerometers  to  measure  the  components  of  the  gravity  gradient  tensor 
from  a  specially  constructed  satellite  with  a  disturbance  compensation  system. 
The  satellite  will  be  launched  to  a  circular,  polar  orbit  of  200  to  250km 


altitude  for  a  six  month’s  mission.  The  accuracy  obtainable  for  the  terrestrial 
gravity  field  is  estimated  to  be  2  to  5  mgal  for  a  resolution  of  150  to  300km 
from  measurement  precision  of  0.01E.  The  expected  launch  data  is  1991  -  1992. 

Simulation  studies  were  carried  out  to  estimate  the  accuracy  and 
resolution  of  the  gravity  field  from  satellite  gradiometry.  Some  of  these 
studies  are:  Breakwell  (1979),  Jekeli  and  Rapp  (1980),  Colombo  and  Kleusberg 
(1983),  Rummel  and  Colombo  (1983),  Rapp  (1985). 

Assuming  a  satellite  gradiometer  with  an  accuracy  of  10“4E  (rms),  Colombo 
and  Kleusberg  (1983)  estimated  the  attainable  accuracies  for  the  gravity  field 
and  for  the  distance  of  the  satellite  to  the  center  of  the  Earth.  A  six  month 
mission  in  polar  orbit  at  200km  altitude  with  data  taken  every  three  seconds, 
would  provide  data  for  computation  of  the  harmonic  coefficients  up  to  degree 
and  order  300  with  less  than  50%  error  and  improve  the  coefficients  up  to 
degree  30  by  up  to  four  orders  of  magnitude  compared  to  existing  values. 

A  simulation  study  suggested  that  an  adjustment  based  on  gradiometer 
data  could  produce  orbital  accuracy  in  radial  distance  10cm  or  better,  if  the 
orbits  are  2000km  high  and  first  an  improved  gravity  model  to  degree  30  is 
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achieved.  This  may  substantially  improve  satellite  altimetry. 

Rummel  and  Colombo  (1983)  in  their  study  assume  a  satellite-gradiometer 
measuring  all  six  second-order  derivatives  of  the  potential.  With  this 
information  it  is  possible  to  separate:  gravity  field  recovery,  altitude  of  the 
satellite,  and  orbit  determination.  This  allows  the  ubo  of  fast  spherical 
harmonic  analysis  like  in  Colombo  (1981a).  The  simulation  study  explains  the 
separation  of  the  orbit  and  the  gravity  parameter  estimation  process.  The 
gravity  field  was  generated  only  from  Zonal  potential  coefficients  up  to  n  = 
300  (The  resulting  gravity  becomes  invariant  in  longitude).  The  actual  orbit 
displacement  and  the  potential  coefficients  were  almost  exactly  recovered  after 


two  iterations. 


Recently  Rapp  (1985)  computed  the  expected  accuracies  for  anomalies  and 
undulations  for  a  six  month  mission  with  a  radial  component  gradiometer 
(*10~*E)  at  an  altitude  of  130km.  The  results  in  terms  of  accuracy  versus 
resolution  are  plotted  on  two  graphs  and  compared  to  an  SST  mission  (*1 
/im/sec.)  and  to  two  models,  OSU  (Rapp)  81  and  GEM  L-2.  On  a  global  average, 
the  SST  mission  improves  the  current  models  by  a  factor  of  10.  The 
gradiometer  mission  improves  the  SST  results  by  about  60%  with  the  exception 
of  the  very  short  wavelengths  (i.e.,  <  20km).  Taken  from  the  anomaly  graph, 
the  gradiometer  could  obtain  an  accuracy  of  *3  mgal  with  50km  resolution,  the 
SST  error  is  three  times  as  large.  From  the  geoid  undulation  graph  for  50km 
resolution,  the  accuracy  of  the  gradiometer  geoid  is  *2cm  versus  the  SST 
geoid  accuracy  of  *18cm  for  the  same  resolution. 


3.4  NASA's  Geopotential  Research  Program  (1982) 


As  it  was  mentioned  in  the  introduction,  the  program  is  divided  into  three 
areas:  a)  interim  field  model  improvements;  b)  geopotential  research  mission; 
and  c)  advanced  mission.  All  three  areas  will  be  discussed  briefly  based 
mainly  on  J.P.  Murphy’s  paper  (1983),  presented  at  the  I.A.G.  Symposium  in 
Hamburg  Germany  (1983). 

The  objective  of  the  "Interim  Field  Model  Improvements"  phase  is  the 
achievement  of  a  "factor  of  two"  improvement  in  the  current  GEM  models.  In 
essence,  this  is  a  rework  and  analysis  of  older  observations  by  utilizing 
current  and  improved  techniques,  updated  force  models,  recently  obtained 
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direct  or  indirect  gravity  observations!  etc.  In  1982  a  Gravity  Field  Workshop 
was  held  at  the  Goddard  Space  Flight  Center  to  specify  the  required  actions 
for  the  achievement  of  the  improvements.  For  the  improvement  of  the  long 
wavelength  features  of  the  models,  the  recommended  steps  were:  1) 

reprocessing  of  selected  laser,  S-band,  Doppler,  and  SST  data  using  improved 
models  for  measurement  corrections  and  for  satellite  motions;  2)  incorporation 


of  orbital  information  obtained  from  the  analysis  of  the  motion  of  synchronous 
and  of  satellites  in  resonant  orbits;  3)  processing  of  additional  Doppler  data 
from  satellites  with  different  orbital  inclinations;  4)  acquiring  and  processing 
laser  data  from  four  satellites  in  a  new  tracking  program;  5)  reprocessing 
available  optical  data  with  improved  star  catalogue  (F5);  and  6)  incorporation 
of  orbit  perturbation  information  on  the  Zonal  coefficients. 

For  the  improvement  of  the  short  wavelength  features,  the 
recommendations  of  the  workshop  were:  1)  improvement  of  data  collection  in 
surface  gravimetry  and  incorporation  of  data  into  the  models;  2)  full  utilization 
of  satellite  altimetry  data  with  state  of  the  art  ephemerides,  ocean  topography, 
tidal  and  other  environmental  corrections.  In  addition,  the  use  of  optimal 
estimation  and  adjustment  techniques  and  of  the  latest  in  computer  technology 
is  recommended. 

The  purpose  of  the  Geopotential  Research  MiBBion  (GRM)  is  to  map  the 
gravity  and  magnetic  fields  with  accuracies  and  resolutions  necessary  for 
requirements  in  geodetic  solid  earth  and  ocean  science  applications. 
Considering  only  gravity  requirements  the  geomagnetic  aspects  of  the  program 
are  omitted  from  this  summary. 

The  GRM  system  will  consist  of  two  satellites  to  be  launched  by  the  space 
shuttle  and  maneuvered  into  a  circular  polar  orbit  of  160km  altitude.  The 
separation  along  the  orbit  between  the  two  satellites  could  vary  between  100 
and  600km.  The  two  satellites  are  equipped  with  a  Doppler  satellite-to-satellite 
tracking  system,  with  a  design  capability  of  one  micron  per  second  relative 
range  rate  measurements.  When  all  the  components  of  the  system  error 
budget  are  taken  into  account,  the  overall  performance  is  better  than  .1 
micrometer  per  second  for  four  second  averages,  one  order  of  magnitude 
smaller  than  the  design  goal  of  one  micrometer  per  second  (Murphy  1983).  Air 
drag  and  other  non-gravitational  forces  are  cancelled  by  the  Disturbance 
Compensation  System  (DISCOS).  this  system  keeps  the  proof  mass  and  the 
satellite  in  a  purely  gravitational  orbit.  The  mission  is  planned  for  six 
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months,  during  which  time  there  will  be  six  tracks  through  each  1*  *  1* 
square.  The  tentative  launch  date  is  1992-1993  (Rapp  1985).  As  it  was  stated 
previously,  the  goal  of  the  mission  is  to  obtain  the  gravity  field  globally  to  *1 
mgal  and  the  geoid  accurate  to  5cm  with  100km  resolutions  for  both.  Most 
significantly,  this  data  set  will  cover  the  currently  inaccessible  and 
gravimetrically  "unknown"  areas.  The  resolution  and  the  accuracy  of  the 
gravity  field  will  allow  the  computation  of  0.5*  *  0.5*  mean  anomalies  globally, 
and  in  turn,  the  expansion  of  harmonic  coefficients,  to  high  degree.  A  better 
gravity  model  will  produce  more  accurate  orbits  for  low  satellites  used  for 
satellite  altimetry  and  for  other  gravity  sensors  as  well.  An  improved  geoid 
will  help  studies  of  sea  surface  topography  and  the  conversion  of  ellipsoidal 
heights,  produced  by  space  techniques,  into  geoid  heights.  In  the  field  of 
geophysics,  the  wavelength  features  of  the  gravity  field,  between  100  and 
1000km  will  help  the  understanding  of  the  mass  variations  in  the  upper  mantle 
(lithospheric  plate  motions).  Shorter  wavelength  gravity  anomaly  features  may 
permit  the  study  of  mantle  convection  and  density  inhoroogeneities  (Taylor  and 
al.  1983). 

The  third  phase  of  NASA’s  Geopotential  Research  Program  Plan  is  the 
"Advanced  Mission".  This  phase  consists  of  advanced  instrument  and  system 
studies  for  the  purpose  of  achieving  the  "ultimate”  in  gravity  and  magnetic 
field  surveys  from  space.  According  to  present  concepts,  the  space  vehicle 
for  these  surveys  will  be  a  Tethered  Satellite  System  (TSS)  and  one  of  the 
candidates  for  the  gravity  sensor  is  the  cryogenic  gradiometer  under 
development  by  Paik  and  al.  (1981).  With  the  potential  of  the  space  shuttle, 
the  tethering  concept  is  under  investigation  for  a  number  of  applications. 
Currently,  NASA,  in  a  cooperative  program  with  Italy,  is  working  on  the 
advanced  development  of  a  tethered  vehicle  and  tether  mechanism  for  possible 
test  flight  in  1987.  A  number  of  potential  applications  are  described  by  Bekey 
(1983).  Gullahorn  and  al.  (1984)  gave  a  presentation  on:  Feasibility  of  Gravity 
Gradient  Measurements  from  a  tethered  Subsatellite  Platform  (EOS  Vol.  65,  No. 
28  Abstract).  NASA’s  objective  with  the  gradient  measurements  from  a 
tethered  satellite  is  the  achievement  of  higher  resolution  in  gravity 
measurements. 


3.S  Data  Processing  and  Adjustment  Techniques 

Most  of  the  algorithms  for  estimation  of  harmonic  coefficients  from  various 
gravity  related  observations  are  based  on  the  standard  least  squares 
adjustment.  From  the  observations  of  each  type  of  measurements!  normal 
equations  are  formed  and  these  are  combined  by  a  weighted  least  squares 
adjustment.  The  increasing  quantity  of  the  observations  and  the  requirements 
for  higher  degree  of  coefficients  made  the  ubo  of  rigorous  least  squares 
procedures  impractical.  The  same  situation  applies  to  the  use  of  the  simple 
least  squares  collocation,  due  to  the  large  matrices  which  must  be  inverted. 

Various  authors  designed  and  used  different  simplified  approaches  (see 
previous  sections  on  global  gravity  models)  to  alleviate  the  computer  burden; 
however,  this  could  not  be  done  without  paying  some  price  for  the  shortening 
of  the  computations.  If  the  expansion  is  truncated,  say  at  degree  36,  which 
can  be  handled  comfortably  by  the  average  computer,  the  finer  detail  of  the 
gravity  field  is  not  adequately  represented  in  the  results.  Therefore,  the 
usefulness  of  the  model  is  limited.  If  the  developments  are  carried  to  higher 
degree  and  order,  say  180,  the  percentage  errors  of  the  coefficients  usually 
attain  100X  around  degree  100  of  the  expansion. 

Recently  several  studies  have  been  accomplished  to  find  solutions  which 
could  improve  the  present  situation  and  together  with  the  improvement  of  the 
computer  technology  (supercomputers)  would  be  capable  to  handle  the 
ever-increasing  amount  of  observations  and  to  produce  accurate  high  degree 
solutions. 

Colombo  (1979,  1981a),  utilizing  the  symmetries  of  data  on  spherical  grids 
and  relations  between  spherical  harmonics  and  Fourier  series,  developed 
efficient  algorithms  for  the  estimation  of  Bpherical  harmonics  to  high  degree. 
Algorithms  are  derived  for  the  evaluation  of  harmonic  coefficients  by  numerical 
quadratures,  and  it  is  shown  that  the  number  of  operations  is  the  order  of  N3 
for  equal  angular  grids  (N  is  the  number  of  lines  of  latitude,  or  the  "Nyquist 
frequency"  of  the  grid).  For  the  error  estimation  of  the  coefficients,  Colombo 
utilizes  the  error  measure  of  least-squares  collocation  and  derives  efficient 
algorithms  for  implementation  on  the  sphere.  The  principle  is  that  for 
"regular"  grids,  the  variance-covariance  matrix  of  the  -  data  consist  of 
"Toepliz-circulant  blocks",  so  it  can  be  both  set  up  and  inverted  very 
efficiently  (Colombo  1981).  Efficient  methods  for  computation  of  covariances 


between  area  means  is  described,  and  numerical  examples  are  demonstrated  in 
the  report. 

A  real  data  set  of  5*  x  5*  mean  gravity  anomalies  obtained  from  Rapp’s 
(1978)  1*  x  1*  mean  values  was  used  for  the  demonstration  of  the  use  of 
optimal  estimators.  The  results  are  compared  with  Rapp’s  coefficients  obtained 
by  numerical  quadratures  up  to  order  36  from  the  original  1*  x  1*  data.  The 
collocation  solution  (Colombo  1981)  follows  very  closely  the  values  of  Rapp’s 
solution,  considered  "true"  values  because  of  the  much  finer  1*  *  1*  grid. 

The  Colombo  (1981)  algorithms  for  harmonic  analysis  have  been 
implemented  for  global  1*  *  1*  anomalies  by  Hajela  (1984).  The  original 
Colombo  algorithms  were  modified  before  their  use  for  the  optimal  estimation  of 
the  coefficients  complete  to  degree  and  order  250  and  for  their  error 
estimates.  This  was  necessary  due  to  the  large  array  size  requirements  for 
the  1*  x  1*  anomalies. 

To  retain  the  efficiency  in  the  computation  of  the  coefficients,  it  was 
required  that  the  variances  of  all  anomalies,  in  a  latitude  band,  be  equated  to 
the  average  value  for  that  band.  This  gives  the  correct  values  only  for  a 

pair  of  coefficients  (Cnni,Snm)  for  any  particular  degree  and  order.  The 

variance  for  each  coefficient  in  the  pair  is  arbitrarily  chosen  equal  (Colombo 
1981,  Hajela  1984). 

The  data  set  used  in  this  study  is  the  set  of  harmonic  coefficients  to 
degree  180  developed  by  Rapp  (1981),  termed  data  set  A.  From  Rapp's 
coefficients,  Hajela  computed  a  global  set  of  1*  *  1*  mean  anomalies.  These 
anomalies  then  were  used  to  derive  potential  coefficient  sets  using  optimal 
estimation  procedures.  For  the  anomaly  data,  several  different  error  estimates 
were  assumed.  Rapp  (1981)  also  computed  a  global  1*  x  1*  anomaly  set  (64800 
values),  adjusting  a  combined  set  of  satellite  derived  coefficients  to  degree  36 
with  a  combined  terrestrial  set  of  1’  *  1*  mean  anomalies.  This  second  set 

was  also  used  in  Hajela’s  study,  termed  anomaly  set  B.  Several  sets  of 

harmonic  coefficients  were  developed  from  the  two  sets  of  global  mean 
anomalies  using  generally  two  sets  of  anomaly  error  estimates.  One  set  termed 
"realistic"  (error  estimate  A)  was  computed  as  the  average  variance  in  each 
latitude  band  from  estimates  of  1*  x  1*  anomaly  errors;  anomalies  computed 
from  potential  coefficients  being  assigned  standard  error  of  30  mgal.  The 
average  anomaly  error  in  a  latitude  band  ranged  from  30  to  4  mgal.  Error 
estimate  B  (idealized  version):  with  no  variation  in  latitude  and  an  assigned 


values  of  *5  mgal.  Several  other  anomaly  error  estimates  were  also  used  in 
various  tests. 

Anomaly  set  B  with  error  estimates  A  and  B  was  used  for  the  computation 
of  potential  coefficients  complete  to  degree  and  order  250.  These  are 
compared  with  Rapp’s  (1981)  set  of  coefficients  (data  set  A).  The  improvement 
of  the  coefficients,  per  degree,  due  to  the  use  of  optimal  estimation  over 
Rapp’s  coefficients  using  de-smoothing  factors  is  about:  8%  at  degree  60,  11% 
at  degree  120,  and  33%  at  degree  180.  In  optimal  estimation,  the  total 

percentage  error  per  degree  does  not  exceed  100%  at  degree  250. 

There  is  no  discontinuity  in  the  degree  variances  of  optimally  estimated 
coefficients  at  degree  60  and  180.  Colombo  (1981a)  tested  his  algorithms  for 

optimal  estimation  of  potential  coefficients  with  a  global  set  of  5*  »  5* 

anomalies.  He  estimated  that  two  hours  of  CPU  times  will  be  required  for  the 
computation  of  coefficients  to  degree  180  with  a  global  set  of  1*  *  1*  anomalies 
on  an  Amdahl  470  V/6-11  computer.  With  a  faster  470  V/8  computer,  Hajela 
used  about  60  minutes  for  the  computation  of  coefficients  to  degree  250. 

S.C.  Bose  and  al.  (1983)  Applied  Science  Analytics,  Inc.,  Canoga  Park, 
California,  examined  the  optimal  estimation  of  harmonic  coefficients  to  high 
degree  from  gravity  anomalies  available  globally.  The  basic  method  is  the 
least-squares  collocation  and  similarly  to  Colombo’s  (1979,1981)  approach, 
consist  of  the  exploitation  of  the  grid  structure  of  the  covariance  matrix, 
thereby  dramatically  reducing  the  computational  requirements.  This  study  is 
an  extension  of  Colombo’s  work.  The  authors  fully  explore  the  structure  and 
effects  of  "rotational  symmetry",  i.e.,  the  sample  need  not  be  uniform  over  all 
latitudes  in  the  grid;  and  consider  meridianal  and  equatorial  symmetries  for 
the  reduction  of  the  computations.  The  study  shows  that  each  parallel  should 
not  have  the  same  number  of  data  samples.  At  high  latitude  the  distances 
between  sample  points  become  less  and  less.  This  crowding  of  sampling  is 
ill-conditioning  of  the  data  covariance  matrix.  When  data  will  be  available 
over  polar  regions,  a  complete  global  solution  of  the  gravity  field  will  be 
necessary;  therefore  the  "thinning”  of  the  sample  data  and  avoiding  the 
ill-conditioning  problem  without  increasing  computational  cost  is  significant. 
The  study  presents  the  concept  and  outlines  the  general  approach  to  the 
solution.  Computational  cost  comparisons  for  different  types  of  symmetries  are 
tabulated. 

Rapp  (1984)  reviewed  some  of  the  methods  that  can  be  used  for  the 


combination  of  satellite  and  terrestrial  gravity  data  for  the  development  of 
high  degree  spherical  harmonic  coefficients.  After  a  short  outline  of  an 
adjustment  process  using  satellite  data  only  and  a  discussion  of  gravity 
anomalies  and  the  boundary  condition,  (defining  gravity  anomalies)  the  report 
outlines  a  combination  method  of  satellite  and  terrestrial  data  by  combining 
the  normal  equations  of  the  two  types  of  data  in  a  least  squares  sense.  In 
this  case  no  corrections  are  applied  to  the  terrestrial  anomalies  due  to  the 
topography;  they  are  interpreted  as  surface  free-air  anomalies  and  the  series 
are  being  evaluated  at  the  surface. 

The  method  for  the  combination  of  satellite  and  terrestrial  measurements  is 
based  on  the  orthogonality  relationship  between  gravity  anomalies  and 
harmonic  coefficients.  The  concept  was  initiated  by  Kaula  (1966)  and  utilized 
by  Rapp  (1978,  1981).  In  the  previous  applications  of  this  method,  two 
assumptions  were  made,  namely  the  spherical  approximation  and  that  the 
gravity  anomalies  are  on  the  surface  of  the  reference  ellipsoid.  Using 
Pellinen’s  (1983)  relationship  between  the  I  and  m  component  of  the  disturbing 
potential  and  gravity  anomaly  on  the  ellipsoid  in  terms  of  the  coefficients  C<m, 
S|m,  Rapp  derived  expressions  for  corrections  to  the  coefficients  computed  by 
the  orthogonality  relationship  from  the  anomalies  on  the  ellipsoid.  Applying 
the  corrections,  the  results  will  be  consistent  with  the  coefficients  obtained 
from  satellite  data. 

The  assumption  that  the  gravity  anomalies  are  on  the  surface  of  the 
ellipsoid  is  also  inaccurate.  They  are  free-air  anomalies  on  the  surface  of  the 
earth.  Therefore,  for  the  application  of  the  procedure  described  in  the 
previous  paragraph,  they  should  downward  continued  to  the  ellipsoid.  To 
avoid  the  problems  regarding  the  validity  of  the  analytical  continuation 
technique,  the  study  suggests  upward  continuation  of  the  surface  anomalies  to 
a  bounding  sphere  (The  Brillouin  sphere)  enclosing  all  the  topography.  This 
would  be  an  alternate  approach  to  both  the  downward  continuation  and  the 
ellipsoidal  problems.  Correction  terms  were  derived  to  the  orthogonality 
formula  for  the  effect  of  the  upward  continuation.  For  the  180  *  180  Rapp 
(1981)  coefficients,  using  a  global  set  of  1*  *  1*  mean  elevations,  correction 
terms  were  computed.  The  comparison  of  the  results  with  the  1981values  gave 
the  percentage  error  in  the  coefficients  caused  by  the  neglect  of  elevation 
effects  (use  of  uncorrected  anomalies)  and  of  the  ellipticity  (use  of  spherical 
approximation  formulas).  From  the  plot  of  the  errors  it  can  be  seen  that  the 


elevation  effect  is  small:  IX  at  low  degrees  and  slowly  rising  to  3X  at  degree 
180.  The  effect  of  spherical  approximation  is  large:  10X  at  degree  75,  20X  at 
degree  130,  and  31X  at  degree  180.  The  anomaly  and  undulation  errors  caused 
by  the  elevation  and  ellipticity  effects  are  tabulated.  The  undulation  error  by 
degree  is  20cm  at  degree  2  and  drops  below  3cm  for  all  degrees  above  20. 
The  cumulative  undulation  error  increases  slowly  from  20  to  34.4cm  at  degree 
180.  Anomaly  errors  gradually  increase  to  3.12  mgal  at  180  degrees.  The 
report  points  out  that  undulation  errors  are  the  largest  in  the  polar  regions 

and  in  high  mountain  areas  or  regions  of  high  vertical  anomaly  gradients 

(Aleutian  region),  when  the  errors  can  reach  80  cm. 

Some  details  and  alternate  adjustment  techniques  are  discussed  for  the 
combination  of  satellite  and  terrestrial  data  using  the  orthogonality 

relationship  and  gravity  anomalies  on  the  bounding  sphere.  The  end  products 
are  the  adjusted  potential  coefficients  and  adjusted  anomalies  on  the  bounding 
sphere.  These  anomalies  can  be  developed  into  high  degree  coefficients  by 
Colombo's  procedure  (Hajela  1984).  The  combination  solution  can  be 
accomplished  with  the  rigorous  formation  of  the  normal  equations.  It  is 

suggested  in  the  report  that  the  large  computer  requirement  can  be  managed 
by  the  use  of  a  supercomputer;  or  Colombo’s  (1981a)  technique,  where  the 
unknown  coefficients  are  ordered  in  a  way  that  the  normal  equations  are  block 
diagonal,  reducing  the  computer  time  substantially.  After  the  combination 
solution,  a  new  set  of  harmonic  coefficients  can  be  computed.  This  expansion 
can  be  developed  to  degree  of  250  (Hajela  1984). 

Another  recommendation  of  the  Rapp  (1984)  report  is  to  use  0.5*  ■  0.5* 
block  size  for  the  mean  anomalies  globally  where  it  is  available.  This  set  of 
anomalies  should  be  upward  continued  to  the  bounding  sphere.  A  global  set 
of  0.5*  *  0.5*  mean  anomalies  can  produce  harmonic  coefficients  to  about  360 
degrees.  Obviously  the  solution  will  be  deficient  in  areas  lacking  data. 
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4.  THE  UTILIZATION  OP  SPHERICAL  HARMONIC  MODELS  IN  GEODESY 


The  availability  of  spherical  harmonic  expansions  to  high  degree  (Rapp 
1978,  1981;  Lerch  and  al.  1981;  Hajela  1984)  made  possible  the  use  of  these 
expansions  for  many  geodetic  and  geophysical  applications  (Tscherning  1983). 
These  expansions,  together  w.th  the  expansion  of  topography  and  of  the 
potential  of  the  isostatically  compensated  topography  (Rapp  1982)  give  many 
computational  advantages  as  compared  to  the  techniques  without  the  use  of 
the  series. 

The  first  obvious  use  of  the  coefficients  of  the  harmonic  expansions  of  the 
gravity  potential  is  the  computation  of  quantities  of  the  gravity  field  such  as: 
gravity  anomalies,  height  anomalies,  components  of  the  deflection  of  the 
vertical,  etc.  For  the  computation  of  these  quantities,  very  efficient 
algorithms  have  been  derived  and  are  in  general  use  today  (Rizos  1979; 
Colombo  1981;  Rapp  1982;  Tscherning  and  Poder  1982;  and  Tscherning  and  al. 
1983).  In  Tscherning  and  al.  (1983),  the  four  above  referenced  programs  are 
briefly  described  and  intercompared  for  results  and  to  obtain  timing 
comparisons. 

The  Rizos  (1979)  program  computed  the  height  and  giavity  anomaly  from 
harmonic  coefficients  for  points  of  a  two  dimensional  evenly  spaced  geographic 
grid.  The  area  of  the  computation  can  be  local,  regional,  or  global. 

The  Colombo  (1981)  computer  program  calculates  the  height  and  gravity 
anomalies  from  potential  coefficients  for  a  global  grid  only  at  specified  latitude 
and  longitude  intervals.  A  subroutine  was  designed  for  the  computation  of 
area  mean  or  point  values.  The  Legendre  functions  are  first  computed  for  the 
grid  interval.  Due  to  the  grid  symmetry  with  respect  to  the  equator,  the 
values  are  computed  for  latitudes  north  of  the  equator.  Then  by  a  Fast 
Fourier  Transform  (FFT),  sums  of  series  are  computed  along  the  latitude  rows. 
This  procedure  is  implemented  by  subroutine  FFTCC  of  the  IMSL  subroutine 
library. 

The  Rapp  (1982)  program  is  for  the  computation  of  the  height  anomaly, 
gravity  anomaly,  gravity  disturbance,  and  the  components  of  the  deflection  of 
the  vertical  from  spherical  harmonic  coefficients.  The  program  has  been 
tested  with  coefficients  up  to  degree  180  and  can  be  extended  higher.  The 
program  calculates  the  above  listed  quantities  on  a  point  to  point  basis. 

The  Tscherning/Goad  Program  (Tscherning  and  al.  1983):  The  original 
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Tscherning  and  Poder  (1982)  program  extends  the  derivatives  of  the  potential 
to  the  second  derivatives  and  uses  the  Clenshaw  summation,  which  has  a 
numerical  advantage  compared  to  the  usual  methods,  consisting  of  a  decrease 
in  the  loss  of  significant  digitB  during  the  summation.  With  this  algorithm  it 
was  possible  to  evaluate  the  sum  of  spherical  harmonic  series  with  degree  and 
order  180  on  a  computer  using  only  10  1/2  significant  digitB  (Tscherning  and 
al.  1983). 

The  original  Algol  programs  given  in  Tscherning  and  Poder  (1982)  was 
translated  into  Fortran  by  Goad,  adding  subroutines  needed  for  the  actual 
computations. 

The  aforementioned  four  programs  have  been  intercompared  by  test  runs 
on  the  Amdahl  470  V/8  computer.  Timing  comparisons  are  tabulated  in 
Tscherning  and  al.  (1983).  Considering  point  to  point  calculation  times,  Rapp's 
and  the  Tscherning-Goad  programs  are  comparable.  Rizos*  is  fastest,  with  0.46 
seconds,  and  Tscherning/Goad  next  with  1.91  seconds  and  the  Rapp  program  a 
distant  third  with  15.59  seconds  (no  provision  is  made  in  the  Rapp  program 
for  data  given  on  a  uniform  longitudinal  grid).  The  time  comparison  of  the 
computation  of  height  anomalies  on  a  global  1*  *  1*  grid  gave  47  seconds  for 
the  Colombo  program  and  66  seconds  for  Rizos.  In  the  other  two  point  by 
point  programs  (Rapp  and  Tscherning/Goad)  the  results  computed  at  one 
latitude  were  extrapolated,  resulting  in  relatively  poor  time. 

It  can  be  seen  that  efficient  algorithms  exist  for  the  computation  of 
gravimetric  quantities  from  spherical  harmonic  expansions  of  high  degree.  For 
limited  or  global  coverage  of  one  or  all  five  quantities,  the  appropriate  best 
fitting  program  can  be  selected. 

Using  high  degree  and  order  coefficients,  the  height  anomalies  can  be 
computed  globally  with  *1.2m  standard  error.  If  the  degree  and  order  of  the 
coefficients  is  increased  to  N  =  360,  and  with  some  improvements  of  the 
coefficients,  the  standard  error  of  the  height  anomaly  could  be  improved  to 
0.5m  (Tscherning  1983).  The  height  anomalies  can  be  used  for  conversion  of 
ellipsoidal  heights  into  normal  heights,  reduction  of  distances,  etc. 

High  degree  spherical  harmonic  expansions  are  used  as  the  reference- base 
for  local  gravity  field  determination,  e.g.  Forsberg  and  Tscherning  (1981), 
Lachapelle  and  Rapp  (1982),  and  Sunkel  (1983).  In  the  future,  the  harmonic 
series  may  replace  the  role  of  normal  potential. 

Due  to  the  easy  computation  of  geodetic  quantities  of  the  gravity  field 
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from  harmonic  series,  the  contribution  of  the  high  degree  expansions  can  be 
subtracted  from  the  point  values  of  these  quantities.  The  residuals  will  have 
empirical  covariance  functions  for  which  the  first  zero  value  is  located  at  a 
spherical  distance  of  about  7*  for  height  anomalies,  45*  for  gravity  anomalies, 
and  20'  for  the  longitude  component  of  the  deflection  of  the  vertical 
(Tscherning,  1983).  Therefore,  the  spherical  caps  used  for  the  evaluation  of 
integral  formulae  will  be  much  smaller,  and  local  collocations  require  a  much 
smaller  data  collection  area.  In  addition,  many  elements  in  the  normal  equation 
matrices  may  be  set  equal  to  zero,  resulting  in  only  about  IX  error  and 
substantial  savings  in  computation. 

Another  beneficial  effect  of  using  high  degree  harmonic  coefficients  is  that 
the  long  wavelength  part  of  the  topographic  effects  are  included  in  the  model. 
For  the  remaining  effects,  the  topography  must  be  considered  only  for  a  short 
distance  from  the  computation  point. 


5.  THE  LOCAL  AND  REGIONAL  GRAVITY  FIELD 


The  description  of  the  gravity  field  of  a  limited  area  (local  or  regional) 
usually  consists  of  the  determination  of  some  functionals  of  the  anomalous 
gravity  potential  in  a  particular  area  and/or  at  selected  points.  These 
functionals  or  gravimetric  quantities  are:  the  components  of  the  deflection  of 
the  vertical  ((, i j);  the  height  of  the  geoid  (N);  the  gravity  anomaly  (Ag)  etc. 
Mathematically  these  functionals  can  be  represented  by  integral  formulae  or  by 
spherical  harmonic  series.  The  two  types  of  formulas  are  equivalent 
theoretically,  but  in  the  practical  work  they  are  different  due  to  the 
differences  in  the  types,  spacing,  distribution,  and  accuracy  of  input  data 
used  for  the  computations  (mean  or  point  anomalies  for  the  integral  formula 
and  geopotential  coefficients  for  the  series  of  spherical  harmonics). 

In  this  section  the  following  topics  are  reviewed:  spectral  properties  of 
data  types;  representation  and  estimation  techniques  of  the  local  gravity  field 
by  integral  formulae  and  collocation;  prediction  of  gravity  outside  the  earth 
from  surface  data  and  improvement  of  local  and  regional  gravity  fields  by 
airborne  gradiometry.  Representative  recent  works  in  each  subject  area  are 
reviewed  and  discussed.  It  is  relied  on  lectures  of  the  Beijing  International 
Summer  School,  China,  1984,  Schwarz,  K.P.  (1984)  and  on  the  published  works 
of  other  authors  in  this  subject  area,  e.g.  Forsberg,  Tscherning,  Svinkel,  Cruz, 
etc. 


5.1  "Data  Types  and  Their  Spectral  Properties" 

In  his  lectures  at  the  Beijing  Summer  School  K.P.  Schwarz  (1984a)  under 
the  above  quoted  title  presents  and  analyzes  a  number  of  topics  relevant  to 
the  representation  of  the  local,  regional  and  global  gravity  fields.  Various 
parts  of  the  lecture  are  summarized  and  freely  quoted  in  the  sequel. 

The  quality  of  local  gravity  field  estimation  depends  on  the  following 
factors: 

-  the  density  of  the  available  data 

-  area  coverage  (global  or  local) 

-  measurement  accuracy 

-  the  sensitivity  of  the  gravimetric  quantity  (functional),  to  be  estimated, 


to  the  given  data  set. 

This  can  be  assessed  by  transforming  the  measurements  into  the  frequency 
domain  and  comparing  the  approximation  obtained  to  the  rigorous 
representation  of  series  of  spherical 

harmonics.  In  view  that  all  functionals  and  data  types  can  be  represented  in 
the  form  of  the  series  of  spherical  harmonics,  the  approximation  error  is  due 
only  to  the  factors  listed  above. 

The  spectral  sensitivity  of  some  specific  functionals  are  evaluated  with 
respect  to  the  following  frequency  ranges: 

2  <  low  <  36 
37  <  medium  <  360 
361  <  high  <  3600 
36001  <  very  high  <  36,000 

The  spectral  sensitivities  of  geoid  height  (N),  gravity  anomaly  (Ag)  and  of  the 
second  order  radial  gradient  of  the  anomalous  potential  (Trr)  are  computed  in 
percentage  of  the  total  value  of  each  spectral  range  (Table  5.1).  The 
computations  were  made  by  the  use  of  the  global  covariance  model  of 
Tscherning  and  Rapp  (1974). 


FREQUENCY  RANGES 

FUNCTIONAL 

low 

■ed. 

high 

v.  high 

N 

99. 2% 

0.8% 

0.0% 

0.0% 

Ag 

22.5% 

41.9% 

32.7% 

2.8% 

Trr 

0.0% 

0.8% 

39.0% 

60.2% 

Table  5.1.  Spectral  sensitivity  of  some  functionals.  After 
K.P.  Schwarz  (1984) 

It  can  be  seen  from  table  5.1  that  for  accurate  determination  of  N  the  low  and 
medium  frequencies  are  required  with  high  accuracy,  and  for  the  Trr  the  high 
and  very  high  frequency  ranges  are  dominant.  There  is  an  exact 
correspondence  between  the  frequency  and  the  spacial  domain;  from  the 
spectral  sensitivity  of  a  functional,  the  spacial  representation  of  the  frequency 
range  in  gridded  form  can  be  determined  (data  density). 
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The  local  gravity  field  estimation  involves  a  specific  data  distribution,  and 
the  question  is  how  well  a  specific  functional  can  be  determined  from  it.  For 
example,  consider  a  set  of  errorless  gridded  data.  The  spectral  representation 
of  this  spacial  data  set  is  a  truncated  series  of  spherical  harmonics.  The 
estimated  coefficients  will  contain  systematic  errors  because  the  data  contains 
the  effects  of  all  frequencies.  This  effect  is  the  aliasing  and  it  is  a  source  of 
errors  of  the  spectral  representation.  If  there  is  only  a  local  data  set,  the 
low  frequencies  cannot  be  properly  determined.  This  is  known  as  spectral 
leakage.  The  integration  formulas  will  have  this  problem  if  the  integration  is 
extended  only  over  a  limited  area.  The  measuring  errors  of  the  data  add 
additional  complications. 

The  steps  to  obtain  a  good  estimation  of  a  local  gravity  field  according  to 
Schwarz’s  lecture  are:  first,  the  spectral  sensitivity  of  the  functional  to  be 
computed  should  be  analyzed;  second,  data  types  which  contribute  most  to  the 
sensitive  part  of  the  spectrum  should  be  selected.  The  aliasing  and  spectral 
leakage  problems  can  be  minimized  by  existing  standard  procedures.  The 
effect  of  measuring  errors  can  be  reduced  if  their  statistical  behavior  is 
known.  In  addition  to  the  methods  listed  above,  the  problems  of  combining 
several  data  types  with  different  spectral  and  statistical  characteristics  muBt 
be  solved.  This  leads  to  the  characterization  of  current  and  future  data 
types. 

The  current  and  future  data  types  are  .characterized  (Table  J3.2)  according 
to:  spectral  resolution,  data  density,  data  coverage,  data  distribution,  and 

noise  spectrum.  These  characteristics  are  useful  for  optimal  data  combination 
for  the  gravity  field  spectrum.  Theoretically,  each  data  type  contains  the 
total  spectrum;  however,  in  practice  the  measuring  process  acts  as  a  bandpass 
filter  limiting  the  range  of  the  spectrum.  Therefore,  a  single  data  type  cannot 
resolve  the  complete  spectrum  and  it  is  necessary  to  combine  different  types 
of  measurements  to  obtain  a  homogenous  spectral  resolution. 

It  can  be  seen  from  Table  5.2  that  the  low  frequency  information  comes 
from  satellite  observations  and  no  other  source  is  in  sight  for  this  frequency 
range.  The  medium  range  is  currently  determined  from  1*  *  1*  mean  gravity 
anomalies  on  land  combined  with  satellite  altimeter  data  over  the  oceans.  This 
frequency  range  will  be  improved  by  satellite  to  satellite  gradiometry,  if  these 
programs  will  become  operational.  The  above  types  are  regular  in  distribution 
and  global  in  coverage,  therefore  their  use  in  integral  formulas  (space  domain) 
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or  in  series  of  spherical  harmonics  (frequency  domain)  should  yield  the  same 
accuracy.  All  the  other  data  types  cover  only  part  of  the  space  domain, 
therefore  integral  formulas  will  give  better  results. 

The  high  frequency  spectrum  is  currently  resolved  by  mean  anomalies  of 
5’  «  5’,  or  deflections  of  the  vertical.  In  the  future,  airborne  gradiometry  and 
inertial  surveys  may  contribute  to  this  frequency  range,  provided  the 
projected  accuracies  of  1  mgal  and  013  will  be  realized.  The  very  high 
frequency  range  currently  is  very  little  known.  Airborne  gradiometry  and 
height  data  may  improve  the  situation  in  the  future.  Some  of  the  very  high 
frequency  range  is  removed  from  gravity  anomalies  by  terrain  corrections 
computed  from  height  data  on  a  grid  (1km  *  1km).  This  removes  the 
frequencies  dependent  on  the  topography.  The  remaining  part  represents  the 
anomaly  variations  caused  by  density  changes.  At  present,  very  little  is 
known  on  the  spectral  power  of  density  changes  in  the  very  high  frequency 
range.  For  the  purpose  of  optimal  spacing  of  gravity  and  height 

measurements,  a  detailed  numerical  analysis  is  recommended  by  K.P.  Schwarz. 
A  spectral  analysis  of  high  density  gravity  and  height  data  would  indicate 
whether  or  not  the  upper  part  of  the  very  high  frequency  spectrum  can  be 
computed  with  sufficient  accuracy  from  height  data  only.  If  yes,  the  interval 
of  gravity  measurements  can  be  larger  than  those  of  height  data,  and  the 
existing  height  data  can  be  used  with  benefit.  Airborne  gradiometer  and 
surface  density  data  would  provide  in  the  future  material  for  Bpectral  analysis 
of  this  frequency  range. 

The  following  parts  of  the  lecture  (Chapters  5  and  6)  discuss  the  concepts 
and  mathematical  treatment  of  spectral  analysis  of  discrete  data  which  apply  to 
local  gravity  field  approximation.  The  one  dimensional  case  is  discussed  first, 
then  the  extension  to  a  two  dimensional  surface.  Correlation  and  spectral 
density  functions  are  described.  The  main  formulas  are  given  with  derivations 
and  some  of  the  literature  is  listed  for  more  details.  For  the  one  dimensional 
case:  Popoulis  (1965),  Bendat  and  Miller  (1982).  For  a  good  discussion  of  all 
essential  formulas  for  the  two-dimensional  planar  case,  I  recommend  Sideris 
(1984),  and  for  the  spherical  case,  Colombo  (1981a).  These  topics  have  also 
been  treated  in  textbooks  and  in  the  related  literature  and  will  not  be 
discussed  here.  It  should  be  noted,  however,  the  use  of  convolutions 
theorems  in  gravity  field  approximations.  The  convolutions  of  two  functions, 
in  general,  can  be  interpreted  as  the  filtering  of  one  function  by  the  other. 
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Examples  are  the  Stokes  and  Vening-Meinesz  integrals.  These  two  integrals 
can  be  considered  as  two  different  filters  applied  to  the  same  gravity  anomaly 
function  (Ag).  The  theorem  states  that  "convolution  in  the  space  domain  can 
be  replaced  by  multiplication  in  the  spectral  domain  and  vice  versa".  All  of 
the  basic  integrals  in  gravity  field  estimation  can  be  expressed  as  convolution 
integrals;  therefore  this  theorem  permits  the  replacement  of  integration  (space 
domain)  by  multiplication  in  the  frequency  domain.  This  technique  is  very 
advantageous  because  recently  gravity  and  height  data  are  available 
frequently  in  gridded  form,  and  Fast  Fourier  Transform  (FFT)  allow 
performance  of  discrete  Fourier  transformations  and  convolutions  much  more 
rapidly  than  integration.  Detail  for  the  Fourier  approach  to  Stokes’  and 
Vening-Meinesz'  integral  are  given  in  Jordan  (1978),  and  details  of  the  terrain 
correction  integral  in  Sideris  (1984) 

The  Spectral  Analysis  of  Global  Gravity  Data  is  discussed  next  in 
connection  with  the  use  of  a  set  of  harmonic  coefficients,  as  a  reference  field 
for  local  or  regional  gravity  data. 

The  concepts  of  spectral  analysis  extended  to  the  sphere,  the  series  of 
spherical  harmonics,  and  their  convergence,  is  widely  covered  in  the 
literature,  e.g.  Moritz  (1980).  The  application  of  the  concepts  to  the  geodetic 
problems  are  given  in  detail  in  Meissl  (1971).  Detailed  numerical  formulas  and 
extension  of  FFT  to  the  sphere  is  discussed  in  Colombo  (1981).  In  principle, 
the  transform  pair  for  the  representation  of  the  anomalous  gravity  field  on 
the  sphere  is  equivalent  to  the  transform  pair  representing  the  gravity  field 
on  the  plane;  however,  the  topology  differences  of  the  two  surfaces  cause  the 
direct  application  of  FFT  techniques  to  the  sphere  to  be  impossible.  In 

Colombo  (1981a),  it  is  shown  that  spherical  harmonics  are  finite  sums  of 
two-dimensional  Fourier  harmonics,  and  this  can  be  used  to  design  efficient 
algorithms  for  spectral  analysis. 

Global  data  available  at  the  present  time  for  spectral  analysis  are: 
Satellite  orbital  perturbations,  satellite  altimeter  data,  and  mean  gravity 
anomaly  data.  The  combination  of  these  three  types  of  very  different 
information  on  the  gravity  field  is  still  not  on  the  desired  optimal  level.  In  a 

typical  procedure,  the  coefficients  are  determined  in  steps,  e.g.  Rapp  (1981). 

The  weighting  is  of  crucial  importance  in  the  combination  of  global  and  local 
gravity  data.  To  obtain  proper  weighting  for  a  global  model,  the 

determination  of  its  error  spectrum  is  necessary;  this  is  difficult  because  the 


global  model  is  a  combination  of  different  types  of  data  with  different  error 
characteristics.  If  a  truncated  series  of  harmonic  coefficients  is  used  as  a 
reference  for  local  approximation  to  obtain  a  proper  weighting  between  the  two 
data  sets,  the  following  questions  arise  (Schwarz  1984): 

(a)  "How  much  of  the  total  spectral  power  of  the  anomalous  gravity  field 
is  contained  in  a  solution  of  degree  and  order  N?" 

(b)  "How  good  is  the  approximation  given  by  a  specific  model;  i.e.,  which 
error  spectrum  is  associated  with  the  global  solution?" 


Practical  procedures  for  weighted  combination  of  harmonic  coefficients  and 
local  data  are  discussed  in  Wenzel  (1982)  and  Sjoberg  (1981). 

Returning  to  the  data  sets  used  in  a  global  solution,  the  characteristics  of 
these  data  is  briefly  discussed.  Orbital  perturbations  can  be  considered  to  be 
upward  continued  geoid  undulations.  Due  to  the  attenuation,  the  short 
wavelength  part  is  smoothed  out  with  the  altitude.  This  effect  and 
measurement  distribution  problems  are  the  principal  reasons  why  only  low 
order  harmonics  can  be  determined  from  this  type  of  data. 

From  altimeter  data,  the  geoid  over  the  oceans  is  obtained,  provided  the 
sea  surface  topography  can  be  neglected.  The  data  are  regular  and  the 
measurement  precision  is  sufficient  for  about  50km  half  wavelength  resolution 
(Cruz  1983).  This  would  be  adequate  for  a  series  expansion  of  360  degree  and 
order,  but  a  similar  set  of  data  is  not  available  over  land  areas.  Currently 
1*  x  1*  mean  anomalies  are  developed  from  the  altimeter  data  and  merged  with 
1*  *  1*  mean  gravity  anomalies  over  land.  Two  particular  problems  emerge 
with  this  data  set:  The  effect  of  averaging  and  aliasing.  Averaging  changes 
the  coefficients.  To  recover  the  spherical  harmonic  coefficients  of  the  point 
gravity  anomaly  function  the  smoothing  by  averaging  has  to  be  reversed. 
This  'desmoothing'  can  be  done  by  optimal  factors  derived  by  Colombo  (1981a). 
The  other  is  the  aliasing,  which  is  not  eliminated  by  "desmoothing".  These 
errors  are  about  50%  of  the  actual  coefficient  at  the  Nyquist  frequency  and  it 
is  larger  above.  Therefore,  there  is  a  need  to  replace  mean  values  with 
smoothed  values  free  of  aliasing.  Simulation  studies  to  used  filters  to 
estimated  mean  values  performed  by  Jekeli  (1981)  eliminated  the  leakage 
problem  for  the  mean  values,  but  not  the  aliasing.  Schwarz's  suggestion  for  a 
possible  way  to  obtain  a  set  of  smoothed  values  free  of  aliasing  is:  "to  use  a 
bandpass  filter  on  the  power  spectrum,  obtain  the  autocovariance  function  by 
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Fourier  Transforming  the  results  to  the  space  domain,  and  estimate  smoothed 
values  in  the  center  of  the  block  using  this  function".  In  addition  to  aliasing, 
missing  mean  anomaly  blocks  and  the  errors  in  the  data  affect  the 
coefficients. 

Recently,  local  data  sets  and  height  data  are  frequently  available  in 
gridded  form,  and  high  degree  geopotential  models  as  reference  fields  are 
available  to  provide  low  frequency  features,  therefore,  the  analysis  of  higher 
frequency  features  in  local  areas  became  possible.  Theoretically  if  a  perfect 
36  «  36  solution  is  available,  a  10*  *  10*  field  of  local  data  centered  on  the 
computation  point  is  required  for  the  resolution  of  all  frequencies  above 
degree  36.  To  improve  the  wavelengths  in  the  medium  range  and  to  help 
spectral  leakage,  a  radius  of  20  degrees  would  be  required  instead  of  the  5 
degree  radius  for  a  10*  *  10*  local  field.  Studies  by  Lachapelle  and  Rapp 
(1982)  indicate,  however,  that  the  extension  of  the  local  area  does  not  improve 
the  long  wavelength  features. 

Here,  like  in  the  case  of  combination  of  global  data  sets,  the  weighting  or 
accounting  for  the  actual  accuracy  of  data  sets  in  a  combination  solution  is  a 
problem.  Wenzel  (1982)  proposed  the  use  of  the  error  covariance  functions  of 
various  data  types  for  the  derivation  of  spectral  weighting  functions  for  an 
optimal  data  combination.  This  approach  is,  in  the  judgment  of  Schwarz,  "most 
promising".  An  advantage  of  the  method  is  that  it  can  be  extended  to  other 
combination  methods,  i.e.  it  is  not  restricted  to  the  integration  method  only. 

Schwarz  analyzed  the  covariance  function  behavior  obtained  by  different 
methods  over  the  same  test  areas,  using  5’  *  5’  mean  anomalies  referenced  to 
GEM-10  5*  area  blocks,  in  North  America.  The  analysis  was  made  by  both  the 
space  domain  and  spectral  methods.  The  difference  in  the  parameters  of 
autocovariance  functions  are  tabulated  for  16  sample  areas  in  Schwarz  (1984). 
Two  dimensional  FFT  methods  have  been  used,  which,  in  addition  to  the 
efficiency,  provide  information  on  the  local  spectrum.  New  results  on  the 
degree  variances  in  the  200  <  n  <2900  are  also  presented. 

5.2  Estimation  of  Gravimetric  Quantities  in  Local  and  Regional  Areas. 

The  estimation  of  the  functionals  of  the  anomalous  gravity  potential  such 
as  deflections  of  the  vertical,  geoid  undulations,  mean  gravity  anomalies,  etc. 


from  the  limited  local  or  regional  information  usually  is  obtained  by  the  use  of 
integral  formulae.  Best  known  examples  are  the  Stokes  and  the 
Vening-Meinesz  integrals.  Theoretically,  the  integration  should  be  extended 
over  the  whole  earth  and  the  gravity  anomaly  should  be  known  at  every 
point.  Practically,  these  conditions  cannot  be  satisfied;  therefore,  these 
integral  formulae  are  usually  modified  to  accommodate  in  combination,  other 
gravity  related  data,  contributing  in  some  form  the  information  outside  of  the 
zone  of  integration.  The  types  of  these  information  are:  harmonic  coefficients 
of  a  series  expansion  of  the  global  field;  information  on  the  topography  in  the 
form  of  digitized  elevation  data,  or  information  on  the  isostatic  compensation  of 
the  topographic  masses.  By  the  use  of  topographic  information  in  the  form  of 
gridded  digital  terrain  models  (DTMs),  the  local  gravity  field  can  be  smoothed 
by  removing  the  effect  of  the  topography  calculated  from  the  DTM.  Terrain 
corrected,  gravity  data  are  of  course  applicable  not  only  in  the  integral 
methods,  but  also  in  the  collocation  technique  for  the  estimation  of  the  gravity 
field. 

Since  spherical  harmonic  expansions  to  degree  and  order  180  are  available 
(Rapp  1978,  1981;  Lerch  and  al  1981)  and  expansions  of  the  topography,  of  the 
rock  equivalent  topography,  and  of  the  isostatically  compensated  topography 
have  also  become  available  (Rapp  1982,  Grasegger  and  Wotruba,  1983),  these, 
together  with  a  fairly  good  coverage  of  DTM  data,  allowed  a  substantial 
number  of  studies  and  numerical  demonstrations  for  the  estimation  of  regional 
and  local  gravity  fields.  Some  of  these  are:  Lachapelle  (1984,  1979); 
Lachapelle,  G.  and  K.P.  Schwarz  (1980);  Lachapelle,  G.  and  A.  Mainville  (1979); 
Tscherning  (1984,  1983,  1983a);  K.P.  Schwarz  (1984a,  1983);  Schwarz  et  al 
(1983);  Schwarz,  K.P.  and  G.  Lachapelle  (1980);  Siinkel  (1984,  1984a,  1983a, 
1983b);  Siinkel,  H.  and  G.  Kraiger  1983);  Forsberg  (1984,  1984a);  Forsberg  and 
Tscherning  (1981);  Rapp  and  Wichiencharoen  (1984);  and  Moritz  (1983,  1980, 
1977). 

5.2.1.  Estimation  by  the  use  of  Modified  Integral  Formulae. 

Lachapelle  (1984),  in  his  lecture  at  the  Beijing  Summer  School,  discussed 
the  modifications  of  the  basic  integral  formulae  of  Stokes  and  Vening  Meinesz, 
and  of  the  topography  integration  formula  for  the  deflection  components,  to 
allow  for  the  combination  of  all  available  data. 

The  free  air  gravity  anomalies  used  in  the  Stokes  and  Vening-Meinesz 


integrals  are  reduced  from  the  earth’s  surface  to  the  geoid  by  the  use  of  the 
gradient  of  the  normal  gravity  (0.3086  mgal/m).  In  mountain  areaB,  the 
vertical  gradient  of  the  actual  gravity  can  differ  as  much  as  0.1  mgal/m.  This 
can  seriously  affect  the  accuracies  of  the  results.  In  the  Molodensky 
approach,  the  measured  gravity  remains  on  the  surface  of  the  earth  and  the 
normal  gravity  is  computed  on  the  telh’roid  by  applying  the  normal  gradient 
upward  from  the  ellipsoid.  The  modified  formulae  (Moritz  1980,  Sec.  48;  and 
Lachapelle  1984)  for  the  height  anomaly  and  for  the  surface  deflection 
components  contain  a  first-order  correction  term  (g,)  to  the  anomaly  on  the 
earth’s  surface  (Ag*  +  g().  It  is  shown  in  Moritz  (1980)  that  the  formula  (g()P 
for  an  arbitrary  point  P  can  be  replaced  by  a  given  terrain  correction 
formula,  if  a  constant  density  of  topography  and  a  linear  correlation  between 
gravity  anomaly  and  height  are  assumed. 

The  terrain  correction  formula  (replacing  the  formula  for  g,)  reveals  that 
the  terrain  correction  is  significant  only  over  a  terrain  with  rough 
topography;  therefore,  the  uncorrected  StokeB  and  Vening  Meinesz  formulae 
(neglecting  the  term  g,)  are  fairly  accurate  approximations  to  the  Molodensky 
approach  in  flat  areas.  The  correction  terms  are  significant  in  mountain 
areas,  particularly  for  the  computation  of  the  deflection  of  the  vertical.  For 
details  see  Moritz  (1980,  section  49).  Numerical  results  for  terrain  corrections 
for  geoid  undulations  are  in  Section  3  of  the  Lachapelle  (1984)  paper.  Four 
different  methods  for  the  combination  of  Stokes’  formula  with  spherical 
harmonic  expansions  are  discussed  and  the  results  of  the  computations  are 
tabulated.  The  spherical  harmonic  models  GEM10B,  and  GEM10C  by  Lerch  and 
al  (1981)  and  two  models  of  Rapp  (1978,  1981)  were  used.  All  except  GEM10B 
are  complete  to  180  *  180.  GEM10B  is  complete  to  order  and  degree  36.  The 
gravity  data  in  form  of  5’  mean  anomaly  coverage  of  North  America  was  used 
(Lachapelle  1978).  Mean  topographic  heights  of  5’  were  used  in  the  Western 
Cordillera  mountain  region  for  the  Computation  of  the  terrain  effect.  The 
undulations  obtained  by  various  methods  are  compared  with  satellite  Doppler 
undulations  (Tables  4,  5,  and  6  of  Lachapelle  1984).  To  show  the  improvement 
in  accuracy  due  to  the  combination  of  surface  gravity  data  with  harmonic 
expansions,  a  comparison  of  Doppler  undulations  and  those  computed  only  from 
geopotential  coefficients  is  given  in  Table  3  of  the  referenced  paper.  The  rms 
values  show  improvements  when  models  complete  to  180  degree  are  used.  Best 
results  of  rms  residuals  of  the  order  of  1  meter  were  obtained  using  Rapp’s 
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1978  model.  The  GEM10B  solution  (36  x  36)  gave  an  rms  agreement  of  1.7 
meters.  Rapp's  1981  180  x  180  set  agreed  to  the  order  of  1.2  meters. 

An  analysis  of  various  modifications  of  Stokes'  function  for  the 
improvement  of  geoid  undulation  computation  accuracy  was  performed  by  Jekeli 
(1981).  Another  recent  study  on  data  requirements  and  accuracy  of  geoid 
undulation  computation  from  gravity  data  is  reported  by  Kearaley  (1984). 

The  deflections  of  the  vertical  are  much  more  sensitive  to  the  gravity 
anomalies  in  the  vicinity  of  the  computation  point  than  undulations. 
Therefore,  deflections  computed  from  harmonic  coefficients  only  will  be  less 
accurate  than  geoid  undulations.  A  comparison  of  deflection  components  of  820 
points  in  Canada  computed  from  harmonic  coefficients  were  compared  with 
astrogeodetic  values.  (Table  7  in  Lachapelle  1984).  It  can  be  seen  from  the 
comparison  of  astrogeodetic  rms  values  with  harmonic  coefficient  derived 
deflection  components  that  the  deflections  from  harmonic  coefficients  only  are 
very  poor.  About  20%  of  the  deflection  signal  comes  from  the  coefficients, 
while  for  the  undulation  signal  the  contribution  of  the  harmonic  coefficients  is 
above  90%. 

Several  methods  for  the  computation  of  the  components  of  the  deflection  of 
the  vertical  from  the  combined  effects  of  gravity  potential  coefficients,  surface 
gravity  and  topographic  data  are  reviewed  in  Lachapelle  (1984). 

The  Vening  Meinesz  formula  was  combined  with  collocation  using  harmonic 
coefficients,  gravity  anomalies,  and  optionally,  astrogeodetic  data  (Lachapelle 
1977).  In  flat  areas  of  Canada,  accuracies  of  1"0  to  1"5  were  obtained  as  rms 
differences  with  astrogeodetic  values.  The  radius  of  the  inner  zone  was  1* 
containing  about  200  5’  x  5*  mean  anomalies.  The  Vening  Meinesz  function  was 
unmodified.  The  modified  function  doeB  not  improve  the  accuracy  in  the  case 
of  deflections,  as  it  is  discussed  in  Jekeli  (1982).  The  computation  software  is 
described  in  Lachapelle  and  Mainville  (1982). 

In  mountain  areas  where  local  topographic  data  are  available,  the  visible 
topography  and  its  isostatic  compensation  can  result  in  about  80%  of  the 
deflection  signal,  provided  that  the  selected  isostatic  model  is  in  fair 
agreement  with  the  reality.  Results  in  the  Swiss  Alps  by  Elmiger  (1969)  and 
in  the  North  American  Western  Cordillera  by  Lachapelle  and  Mainville  (1979) 
confirm  the  estimate. 


In  mountain  areas,  where  astrogeodetic  deflections  exist  in  close  intervals, 
a  combination  of  these  deflections  with  topographic-isostatic  ones  improves  the 
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accuracy.  Elmiger  (1969)  obtained  about  1"5  to  Z!0  in  the  Alps. 
Topographic-isostatic  deflections  alone  may  be  significantly  biased.  Lachapelle 
(1975)  combined  topographic-isostatic  deflections  with  global  data,  such  as 
harmonic  coefficients  of  the  geopotential  for  the  purpose  of  reducing  large 
scale  systematic  effects.  The  deflection  components  were  the  sums  of  the 
contributions  of:  the  low  degree  potential  coefficients,  the  corresponding 

coefficients  of  the  isostatic  reduction  potential,  and  the  topographic-isostatic 
deflection  components  computed  from  the  detailed  topography  using  the 
topography  integration  formula.  Tables  (9  and  10)  in  Lachapelle  (1984), 
extracted  from  Lachapelle  and  Mainville  (1979)  give  some  results  obtained  by 
this  method  in  the  Canadian  and  New  Mexico  (White  Sands)  mountain  areas  at 
selected  astrogeodetic  deflection  stations.  The  rms  differences  with 
astrogeodetic  values  are  given  for  the  topographic-isostatic  components  and 
for  the  total  components.  In  the  Canadian  Cordillera,  mean  differences 
between  the  topographic-isostatic  and  astrogeodetic  components  are  -0128  (() 
and  01 21  (ij)  and  between  the  total  and  astrogeodetic  components  are  0126  (() 
and  -0122  (?/).  This  indicates  that  the  harmonic  coefficients  of  the  geopotential 
and  isostatic  reduction  potential  do  not  remove  significantly  regional  trends. 
In  the  White  Sands  area,  the  predicted  total  deflection  components  agree 
better  with  the  astrogeodetic  values  than  with  the  topographic-isostatic 
components.  Thus,  the  geopotential  and  the  isostatic  reduction  potential 
coefficients  contribute  substantially  to  the  removal  of  regional  biases.  The 
mean  differences  between  the  topographic-isostatic  and  astrogeodetic 
components  are:  -1"37  (()  and  -1"92  (7j).  These  are  reduced  to  (7151  ({)  and 
-01 02  (■»/)  between  the  total  combined  and  astrogeodetic  deflection 
components. 

5.2.2.  Comparison  of  Gravimetric  and  Satellite  Derived  Undulations 

Gravimetric  geoid  undulations  were  computed  by  three  different  methods 
and  compared  with  geoid  undulations  obtained  from  satellite  Doppler  data  at 
the  same  sites  (Rapp  and  Wichiencharoen  1984).  Twenty  Doppler  sites  were 
selected  from  the  65  in  the  United  States  from  the  Lachapelle  (1979)  study,  for 
the  comparison  of  Doppler-derived  and  gravimetric  geoid  undulations  in  North 
America.  Ten  of  these  sites  were  in  the  mountainous  western  United  States, 
where  large  residual  differences  were  found  between  the  Doppler  derived  and 
gravimetric  undulations.  The  other  10  stations  were  in  the  relatively  flat 
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eastern  part  of  the  United  States.  The  original  Doppler  undulations  were 
corrected  for  a  -0.4  ppm  scale  change  and  for  a  4m  Z  axis  bias.  The 
gravimetric  geoid  undulations  were  computed  by:  a)  combining  Rapp  (1981) 
potential  coefficients  to  degree  36  with  uncorrected  1*  *  1*  mean  free  air 
gravity  anomalies  in  a  10*  cap  surrounding  each  station,  b)  using  Helmert’s 
second  condensation  procedure  for  anomaly  reduction  where  terrain 
corrections  and  topographic  indirect  effects  are  computed  (10  stations  in  the 
western  U.S.),  c)  the  use  of  potential  coefficients  to  degree  180  (Rapp  1981). 
The  components  of  the  gravimetric  undulations  and  the  differences  between 
these  and  the  Doppler-derived  values  are  tabulated.  The  effect  of  the  terrain 
correction  is  of  the  order  of  1.8m  and  of  the  indirect  effect  is  15cm.  The 
mean  difference  between  Doppler  and  gravimetric  undulations  for  the  10 
western  stations  excluding  the  terrain  correction  and  indirect  effect  was  1.65m 
with  standard  deviation  of  the  difference  of  *0.96m.  With  the  corrections,  the 
mean  systematic  difference  decreased  to  0.1m  *  0.93m. 

In  the  eastern  U.S.  stations,  where  no  terrain  corrections  are  applied,  the 
mean  difference  between  Doppler  and  gravimetric  undulations  is  0.42  *  0.55m. 

The  spherical  harmonic  coefficients  to  degree  180  gave  results  close  to 
those  when  uncorrected  1*  *  1*  mean  anomalies  were  used.  For  the  10 
stations  in  the  western  U.S.,  a  difference  of  1.27  *  1.07m  was  obtained.  For 
the  eastern  U.S.  stations  it  was  0.45  *  0.62m. 

5.2.3.  On  the  Accuracy  of  Height  Anomalies  and  Deflections  of  the  Vertical 


Obtained  from  Combination  of  Mean  Anomalies  and  Harmonic 


Coefficients. 

The  accuracies  of  height  anomalies  and  deflections  of  the  vertical  obtained 
by  the  combination  of  harmonic  coefficients  and  mean  gravity  anomalies  was 
investigated  by  Heck  (1983). 

It  is  known  that  the  applications  of  the  Stokes  and  Vening  Meinesz 
integrals  require  continuous  gravity  data  over  the  entire  earth.  Since  the 
existing  global  coverage  of  gravity  anomalies  is  neither  continuous  nor 
homogeneous,  in  practice  the  integration  is  limited  to  a  spherical  cap  around 
the  computation  point.  Because  the  gravity  data  usually  is  known  in  terms  of 
mean  anomalies  of  surface  blocks,  the  integrals  are  replaced  by  summations. 

The  outer  zones  are  replaced  or  represented  by  a  geopotential  model. 
Four  types  of  errors  in  such  combinations  are  considered  in  the  Heck  (1983) 
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paper:  1)  error  due  to  the  lack  of  higher  degree  coefficients,  2)  errors  in 

the  used  set  of  coefficients,  3)  error  due  to  the  loss  of  information  by  using 
finite  block  size,  and  4)  errors  of  the  mean  anomalies.  Error  analyses  reveal 
that  at  the  zeros  of  the  kernel  functions  within  the  spherical  integrals,  the 
error  functions  show  local  minima.  This  was  the  motivation  for  modifications 
of  the  kernel  functions  by:  Molodensky  and  al  (1962),  Wong  and  Gore  (1969), 
Meissl  (1971),  and  Colombo  (1977).  The  various  approaches  have  also  been 
studied  by  Jekeli  (1981),  Rapp  (1980),  Chen  (1981),  Fell  and  KaraBka  (1981)  and 
Heck  and  Gruninger  (1983). 

The  investigations  of  Heck  show  that  the  errors  caused  by  the  first  two 
error  sources  listed  above  show  distinct  minima  if  the  integration  radius  is 
confined  to  the  zeros  of  the  kernel  functions.  This  characteristic  suggests 
that  the  Wong  and  Gore  modification  procedure  be  used  for  combining  gravity 
anomalies  of  a  spherical  cap  and  geopotential  coeff'  ients.  It  is  emphasized 
that  the  integration  radius  be  extended  exactly  to  the  first  zero  of  the 
modified  Stokes  and  Vening  Meinesz  kernel  functions. 

Using  a  recent  geopotential  model  and  mean  anomalies  of  2  to  10km  side 
length,  absolute  undulations  have  a  global  rms  of  *1  to  3cm;  relative 
undulations  between  points  50  to  100km  apart  have  an  accuracy  of  a  few 
centimeters.  The  deflections  of  the  vertical  have  an  accuracy  to  an  order  of 
magnitude  of  012  to  013  with  2*2  km2  mean  anomaly  elements.  The  above 
numerical  estimates  have  been  obtained  by  the  use  of  hypothetical  error 
models;  therefore,  the  results  must  be  considered  to  represent  only  the  order 
of  magnitude  of  the  errors. 

The  study  emphasized  the  importance  of  the  consideration  of  terrain 
corrections  for  gravity  anomalies  in  rugged  mountainous  terrain.  Substantial 
systematic  errors  may  occur  using  uncorrected  anomalies  for  computation  of 
gravimetric  quantities. 


5.2.4.  Estimation  of  Gravity  Field  by  Collocation. 

The  integral  procedures,  such  as  Stokes’  and  Vening  Meinesz’,  use  one 

type  of  data  for  the  approximation  of  other  functions  of  the  gravity  field. 

Different  types  of  data  are  frequently  available,  containing  useful  information 
regarding  the  gravity  field.  As  it  is  well  known,  the  collocation  method  is 

capable  of  using  gravity  dependent  heterogeneous  data  to  predict  any  other 

gravity  field  quantity.  It  is  also  known  that  there  is  a  strong  correlation  of 
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the  free  air  anomalies  with  elevation.  This  correlation  results  in  a  trend  in 
the  deflections  of  the  vertical,  in  the  gravity  anomalies,  etc.  which  must  be 
removed  before  the  use  of  the  data  for  interpolation  or  before  applying  them 
in  collocation. 


5.2.4. 1.  Forsberg  and  Tscherning  (1981).  This  paper  describes  various 
methods  for  the  computation  of  terrain  effects  applicable  to  the  estimation  of 
gravity  field  quantities  by  collocation.  The  various  "reduction"  methods  were 
tested  in  two  mountainous  1*  x  1*  areas  in  New  Mexico.  Both  areas  are 
characterized  by  a  North-South  mountain  chain  800-1500m  above  the  plateau 
which  is  about  1200-1400m  high.  In  these  areas,  known  gravity  anomalies  and 
deflections  of  the  vertical  were  predicted.  The  predictions  were  made  by  the 
use  of  stepwise  collocations,  described  by  Tscherning  (1974,  1978)  and 
previously  used  in  Tscherning  and  Forsberg  (1978).  In  the  first  step,  the 
GEM10B  set  of  harmonic  coefficients  was  used.  In  the  second  step  a  set  of 
1*  x  1*  mean  gravity  anomalies  wbb  used.  In  the  third  step,  a  local 
approximation  was  computed  using  anomaly  spacing  of  about  10km  and  a  few 
deflections  40km  apart.  For  this  step,  local  empirical  covariance  functions 
were  estimated  for  anomalies  and  deflections. 

The  predictions  were  carried  out  with  the  unchanged  original  data,  and 
with  the  data  reduced  with  the  following  methods: 

a)  fixed  sector  topographic/isostatic  reduction  for  a  6*  x  7*  area  around  the 
test  areas,  Airy  isostatic  model,  crustal  thickness  32km  and  density 
contrast  (crust/mantle)  0.4  g/cro3 

fixed  sector  (6*  x  7*)  reduction  for  the  visible  topography  above  sea  level 
residual  terrain  model  (RTM)  reduction  with  mean  topography  defined  by  a 
bilinear  interpolation  in  a  30’  x  30’  mean  height  grid,  residual  topography 
taken  into  account  in  a  fixed  3*  x  4*  sector. 

RTM  reduction  with  15’  x  15’  mean  height  grid,  residual  topography  taken 
into  account  out  to  a  distance  of  only  60km. 

Both  the  point  observations  and  mean  gravity  values  were  terrain  corrected. 
The  RTM  effects  on  the  mean  anomalies  were  very  small.  comparing  the 
original  and  the  reduced  observations,  it  is  apparent  that  the  major  par*  of 
the  change  in  the  observed  quantities  is  due  to  the  effect  of  the  isostatically 
compensated  topography.  The  choice  of  the  isostatic  parameters  have  little 
effect  on  the  variation.  It  is  pointed  out  that  some  isostatic  effect  should 
always  be  included  to  avoid  the  bias  on  gravity  anomalies  (and  height 
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anomalies)  which  occur  when  only  the  topography  is  removed.  (Refer  to  Table 
1  of  Forsberg  and  Tscherning  (1981).  A  total  of  110  deflection  pairs  and  150 
gravity  free  air  anomalies  were  predicted  and  compared  to  the  known 
(observed)  values.  The  mean  values  of  the  predicted  deflection  components 
and  gravity  anomalies  in  the  two  1*  *  1*  areas  are  tabulated  with  their 
standard  deviations.  The  values  are  tabulated  in  three  variations  according  to 
the  point  observations  used:  A)  gravity  plus  deflections,  B)  gravity  alone, 
and  C)  no  point  values,  only  the  1*  x  1*  mean  anomalies  and  GEM10B 
coefficients  were  used.  The  tabulation  demonstrates  the  improvement  of  the 
prediction  when  the  effect  of  the  terrain  is  considered.  The  rms  error 
decreased  with  a  factor  of  close  to  3  when  terrain  corrected  data  is  used.  It 
is  also  shown  that  better  deflection  results  are  obtained  using 
topographic/isostatic  reduction  without  any  local  gravity  data  (variation  "C" 
above)  than  by  the  use  of  all  the  gravity  data  without  terrain  effect 
correction.  It  is  concluded  that  it  is  possible  to  predict  deflections  of  the 
vertical  and  gravity  anomalies  in  areas  of  rugged  terrain  with  accuracies  of  1" 
and  3-4  mgal  respectively  from  anomaly  data  spaced  6’  apart,  and  when  the 
gravity  field  is  smoothed  by  terrain  corrections  computed  on  the  basis  of 
0.5  x  0.5  minute  point  heights. 

If  only  rough  height  information  is  available  (such  as  5’  x  5’  mean 
elevations),  substantial  improvements  occur  in  the  predictions.  An  example  is 
the  geoid  prediction  in  Greenland  by  Forsberg  and  Madson  (1981).  The  study 
shows  that  the  collocation  method  yields  results  in  areas  of  rough  topography 
comparable  in  accuracy  to  the  results  obtained  in  flat  areas. 

5. 2. 4. 2.  "The  Geoid  of  Austria"  by  Sunkel.  Another  example  of  the  use  of 
global  geopotential  and  digital  terrain  models  in  combination  with  local  gravity 
information  for  the  computation  of  the  regional  geoid  is  "The  Geoid  of  Austria” 
(Sunkel  1983). 

For  this  work,  the  following  data  was  used:  The  Rapp  180  x  180 

geopotential  model,  1*  x  1*  mean  values  of  the  global  Digital  Terrain  Model 
(DTM),  a  20"  *  20"  digital  terrain  model  of  point  values  interpolated  from  a 
1:500,000  topographic  map  of  Austria  [a  more  sophisticated  program  for  digital 
models  of  topography  and  rock  densities  is  described  in  Steinhauser  and  al 
(1983).),  and  521  observed  astrogeodetic,  deflections  of  the  vertical.  The 
deflections  are  the  results  of  an  extensive  program  started  in  1975  with  new 
instruments.  The  formal  accuracies  of  *  0.2  -  0.3  seconds  were  changed  to 


*  0.5  -  0.7  seconds  as  a  priori  estimates  in  the  collocation.  The  deflections, 
originally  on  a  local  system,  have  been  converted  to  the  geocentric  Doppler 
system  of  the  Naval  Weapons  Lab  (NWL  9D)  and  to  the  ellipsoid  of  the  Geodetic 
Reference  System  1980. 

The  harmonic  coefficients  of  the  gravitational  potential  of  the  topography 
and  its  isostatic  compensation  were  determined  from  the  DTM  data.  Then  the 
coefficients  of  the  Rapp  geopotential  model  and  the  topographic-isostatic 
coefficients  have  been  subtracted  from  each  other  and  the  effect  on  the 
surface  deflections  determined. 

The  topographic-isostatic  data  reduction  was  calculated  in  four  steps.  In 
the  inner  zone  (T  <  2km),  the  exact  parallelpiped  formula  was  used,  and  in  the 
second  zone  (2  <  r  <  38km), a  much  simpler  formula.  In  zones  3  and  4,  the 
mass  point  formula  replaced  the  exact  method.  The  effect  of  the  remote  zones 
(38  <  T  <  150km)  and  (T  >  150km)  were  determined  by  interpolation  from  a 
grid  of  15’  *  15’  and  45’  *  60’  respectively. 

The  two  parameter  "attenuated  white  noise  model"  of  Jordan  and  Heller 
(1978)  was  chose  for  the  model  of  covariance  functions.  Height  anomaly 
covariance  functions  have  been  determined  for  5  regions  within  the  whole 
area.  A  variance  of  0.05m*  and  a  correlation  length  of  35  km  was  obtained  for 
the  area  of  investigation. 

For  the  collocation  process,  a  stepwise  collocation  solution  was  chosen, 
with  the  following  arrangements: 

"At  step  zero,  a  subset  of  30  approximately  homogeneously  distributed 
vertical  deflections  have  been  processed.  For  this  subset,  the  covariance 
matrix,  its  inverse,  and  the  solution  vector  have  been  calculated.  With  this 
data, the  non-processed  vertical  deflections  have  been  predicted  and  compared 
to  the  measurements.  Those  (10)  vertical  deflections,  which  showed  the 
largest  discrepancies,  were  used  as  additional  data  in  the  next  step.  In  thiB 
way,  the  covariance  matrix  and  its  inverse  and  the  solution  vector 
successively  improved  using  block  partitioning  methods." 

It  is  interesting  to  note  the  improvement  of  the  prediction  by  the  stepwise 
approach  with  the  increase  of  the  amount  of  data.  The  rms  vertical  deflection 
error  is  almost  constant  (about  210)  for  up  to  25%  of  the  total  data,  with  50% 
of  the  data  the  error  drops  to  1"0,  and  with  90%  the  error  is  about  013,  It  is 
also  noted  that  the  external  prediction  error  estimates  agree  with  the 
collocation  prediction  error  estimates  for  almost  all  points  with  20%.  This  is 


an  indication  of  the  reliability  of  collocation  error  estimates. 

The  collocation  contribution  to  the  geoid  height  or  the  height  anomaly 
have  been  predicted  for  a  (4,X)  grid  of  3'  in  latitude  and  5’  in  longitude  (2240 
points).  Due  to  the  long  CPU  time  for  the  estimation  of  the  prediction  error, 
this  was  computed  only  for  a  few  sample  points.  The  estimated  rras  errors  of 
the  relative  geoid  heights  range  between  *4  and  *14cm.  The  representative 
estimate  for  the  whole  area  was  *8cm. 

By  adding  the  results  obtained  from  the  topographic-isostatic  reduction 
and  from  the  geopotential  model  to  the  contribution  of  the  collocation,  and 
applying  the  indirect  effect,  the  geoid  height  and  the  actual  height  anomalies 
are  obtained.  The  contributions  from  the  topographic-isostatic  model  and  from 
the  geopotential  model  vary  from  41.5  to  47.5m,  and  the  collocation  contribution 
is  between  *0.5m  after  removal  of  a  trend  on  the  order  of  3m.  The  effect  of 
the  analytical  continuation  on  the  isostatically  and  gravity  model  reduced 
potential  is  of  13cm  in  the  Central  Alps.  The  maximum  difference  between 
height  anomalies  and  geoidal  heights  is  located  at  the  "Grosslockner”  massif 
(3800m),  and  has  a  value  of  35cm,  which  is  in  excellent  agreement  with  the 
formula  given  in  Heiskanen  and  Moritz  (1967,  p.325). 

A  number  of  additional  papers  were  presented  on  the  geoid  in  Austria  at 
the  XVIII  General  Assembly  of  IUGG/IAG  in  Hamburg,  August  1983: 
Bretterbauer  (1983),  Steinhauser  and  al  (1983),  and  Erker  (1983).  A 
comprehensive  description  of  the  work  on  "Geoid  in  Austria"  is  given  in 
German  by  Lichtenegger  and  al  (1983), 

Moritz  (1983)  in  his  report  on  "Local  Geoid  Determination  in  Mountain 
Regions"  reviews  the  methods  for  regional  determination  of  the  geoid  or  of 
height  anomalies  from  deflections  of  the  vertical  or  by  combination  of 
deflections  and  anomalies  by  collocation.  Basic  definitions  and  geometry  are 
reviewed,  geoid  determinations,  reduction  for  curvature  of  the  plumb  line  and 
topographic-isostatic  reduction  of  vertical  deflections  are  discussed  in  the 
classical  and  modern  (Molodensky)  approaches.  The  application  of  the 
collocation  is  outlined  and  a  summary  of  the  geoid  computation  for  Austria 
(SCinkel  1983)  is  given. 

At  the  IAG  Symposium  in  Hamburg  during  the  XVIII  General  Assembly  of 
IUGG/IAG,  Aug.  1983,  a  paper  was  presented  by  Gurtner  and  Elmiger  (1983) 
summarizing  the  work  done  in  Switzerland,  in  the  past  15  years  and  currently, 
regarding  geoid  and  vertical  deflection  computations. 


5. 2. 4. 3.  Other  Gravity  Field  Estimations  by  Collocation.  An  algorithm  for  the 
prediction  of  free  air  anomaly  is  presented  in  Sunkel  and  Kraiger  (1983).  The 

algorithm  was  tested  with  4  data  sets  in  Austria.  The  computation  method 

used  was  least-squares  collocation  with  parameters  assuming  a  linear 
correlation  between  terrain  corrected  free-air  anomaly  and  elevation.  The 
computations  showed  that:  a)  point  free  air  anomalies  can  be  predicted  with  a 
standard  deviation  of  <3  mgal  (or  better)  from  terrain  corrected  free-air 

anomalies  with  data  density  of  better  than  5  data  per  100  km2,  b)  the 

collocation  determined  Bouger  factor  deviates  from  the  standard  value  of  0.112 
mgal/m  less  than  5%  with  an  rms  error  of  *3%;  c)  the  covariance  function  of 
terrain  corrected  and  trend  reduced  free-air  anomalies  agrees  very  well  with 
the  gravity  anomaly  covariance  function  obtained  from  topographic-isostatically 
reduced  vertical  deflections. 

The  actual  accuracies  of  the  predictions,  by  comparing  predicted  values  to 
non-process  measurements  were:  in  the  north-east  foothills  of  the  Alps  the 
prediction  error  for  the  free-air  anomalies  was  ‘3-4  mgal;  in  the  eastern 
foothills  ‘2  mgal;  and  in  the  central  Alps  ‘10  mgal.  The  ‘10  mgal  error  is 
probably  due  to  the  incomplete  terrain  correction  of  free  air  anomalies. 

Tscherning  (1984)  reviews  the  theory  and  the  steps  of  the  implementation 
of  the  use  of  least  squares  collocation  for  computation  of  the  disturbing 
potential  and  related  parameters.  It  is  shown  how  to  reduce  the  number  of 
observations  under  certain  conditions  in  local  solutions.  The  effect  of  the 
smoothing  of  the  gravity  field,  achieved  by  subtracting  out  the  contribution  of 
the  topography  and  eventually  of  the  geological  structures,  is  one  of  the 
topics  discussed  in  detail.  The  computer  processing  is  described,  dividing  it 
into  separate  steps.  Each  step  is  illustrated  in  Tscherning  (1982).  It  is 
shown  that  the  collocation  method  is  a  very  powerful  tool  to  solve  geodetic 
problems  and  has  two  major  advantages,  one  being  that  it  can  process 
heterogeneous  data,  and  the  other  being  its  ability  to  compute  error  estimates. 
Tscherning  (1981)  compared  collocation  with  other  methods  and  concluded  that 
other  methods  "having  a  sufficiently  solid  theoretical  basis"  would  give 
comparable  results. 

Forsberg  (1984)  treats  and  analyzes  a  number  of  topics,  all  related  to  the 
modeling  of  the  local  and  regional  gravity  field.  The  first  topic  is  the 
combined  collocation  and  geophysical  inversion  procedure  for  modeling.  The 


term  "gravity  field  modeling"  in  contrast  to  the  geodetic  interpretation  of 


representing  the  external  potential  of  the  earth,  in  geophysics  "stands  for  the 
process  of  determining  internal  density  distribution  of  the  earth  consistent 
with  the  observed  outer  field".  The  geophysical  modeling  of  the  gravity  field, 
therefore  turns  out  to  be  the  determination  of  the  external  gravity  potential 
due  to  density  distributions.  Our  current  knowledge  of  density  anomalies 
contributes  only  to  the  shorter  wavelength  part  of  the  gravity  field.  The 
longer  wavelength  parts  can  be  obtained  from  spherical  harmonic  expansions, 
termed  "external"  modeling  in  the  report.  The  "external"  and  "internal" 
modeling  could  be  carried  out  simultaneously  by  combined  versions  of 
collocation  and  geophysical  inversion.  The  results:  the  external  gravity  field 
and  density  information  could  be  obtained  simultaneously.  The  theoretical 
background  is  outlined  regarding  density  anomalies  and  an  outline  is  given  in 
principle  on  the  combination  of  inversion  and  collocation  for  gravity  modeling 
to  include  density  (non-gravity)  information  in  gravity  field  approximation. 

The  basic  terrain  reduction  methods  are  reviewed  in  the  follow-on 
sections.  It  is  pointed  out  that  the  conventional  gravimetric  "terrain 
correction"  is  not  a  terrain  reduction;  the  term  is  used  in  the  report  for  a 
small  non-linear  correction  to  the  Bouger  reduction.  A  FORTRAN  program  for 
the  computation  of  terrain  reductions  for  gravity  anomalies,  geoid  undulations, 
and  deflections  of  the  vertical  is  listed  in  the  appendix. 

The  accuracy  of  the  "linear  approximation"  (an  approximative  terrain 
reduction),  used  for  error  studies  and  FFT  techniques,  is  investigated  and  it 
is  found  that  the  approximation  is  usually  acceptable  for  theoretical  models 
and  for  actual  data.  By  assuming  the  validity  of  linear  approximation  FFT 
methods  for  terrain  effect  computation  are  outlined. 

Error  studies  are  made  on  the  Digital  Terrain  Model  (DTM)  resolution 
requirements.  Formulae  and  error  diagrams  are  derived  for  the  computation 
of  rms  errors  of  terrain  effects  on  gravity  anomalies,  deflections  of  the 
vertical,  and  height  anomalies  on  the  basis  of  spacing  of  elevation  data  and  on 
the  covariance  functions  of  the  topography.  The  effect  of  topography  is 
studied  by  investigation  of  empirical  covariance  functions  in  five  different 
areas  of  the  United  states.  These  studies  resulted  in  actual  empirical 
information  on  covariance  functions  for  topography  and  gravity.  In  addition, 


information  on  magnitude  of  terrain  corrections  and  of  RTM  geoid  effects  are 
obtained.  Information  is  given  on  the  degree  of  smoothing  obtained  from 
terrain  reductions  of  actual  gravity  data. 
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Finally,  from  isostatic  reductions  of  satellite  altimetry  data  in  two  areas 
of  the  Pacific,  a  relationship  is  derived  between  ocean  geoids  by  altimetry  and 
those  computed  from  bathymetric  data. 

5.2.5.  Prediction  of  Gravity  Outside  the  Earth  From  Surface  Data. 

The  prediction  of  the  gravity  disturbance  vector  at  high  altitudes  is 
investigated  for  several  technical  reasons,  such  as  high  altitude  gravity 
experiments,  inertial  navigation,  etc.  The  task  is  to  estimate  the  gravity 
disturbance  vector  components,  gravity  anomalies,  or  other  functions  of  the 
gravity  potential  at  altitudes  from  surface  gravity  data,  with  a  theoretically 
valid  and  computationally  manageable  and  feasible  procedure.  In  addition,  it 
is  desired  to  account  for  all  error  sources  involved  during  the  process, 
determine  the  propagation  of  these  errors  into  the  results,  and  obtain  at  the 
end  a  reliable  error  estimate  for  the  final  products.  Basic  steps  of  possible 
procedures  such  as  gravity  reduction  methods,  computations  of  mean 

anomalies,  an  techniques  for  upward  continuation  of  gravity  are  well  known, 
however,  the  actual  work  performed  in  this  area  to  date,  is  experimental  in 
nature. 

Some  of  the  recent  studies  on  this  subject,  selected  for  more  detailed 
description  are  Sunkel  (1984a,  1981),  and  Cruz  and  Laskowsky  (1984). 

5. 2. 5.1.  Prediction  of  the  Gravity  Disturbance  at  Altitudes.  Sunkel  (1984a) 
examines  four  approaches:  The  Green’s  approach,  by  representing  the 

external  disturbing  potential  with  a  linear  combination  of  a  single  and  double 
layer  potential;  the  surface  layer  approach;  a  combination  of  analytical 
continuation  down  to  sea  level  with  the  Pizzetti-Stokes  method;  and  the 

collocation  approach.  Each  method  has  its  advantages  and  disadvantages.  The 
integral  formulae  have  the  common  assumptions  that:  a)  gravity  anomalies  are 
given  at  every  point  of  the  earth,  b)  only  gravity  anomalies  are  used  as  data, 
and  c)  all  data  are  exact.  Obviously,  these  assumptions  are  never  met.  The 
advantage  of  the  integral  formulae  is  that  the  inversion  process  is  performed 
analytically,  resulting  in  a  very  fast  and  simple  process  for  the  computation  of 
gravity  field  quantities.  Conceptually,  the  Least  Squares  Collocation  (LSC)  is 

superior  to  all  other  methods.  It  combines  heterogeneous  input  data, 

downward  continuation  of  surface  data,  smooth  interpolation,  upward 


continuation,  noise  filtering,  and  estimation  of  parameters.  "Theoretically  the 
whole  apparatus  of  LSC  can  be  directly  applied  to  the  unreduced  data  as 
input  yielding  any  estimable  gravity  field  quantity  as  output"  (Sunkel  1984a). 
Theoretically  LSC  seems  to  be  the  ideal  tool  for  the  tasks  under  discussion, 
however,  the  solution  of  a  very  large  system  of  equations  with  a  fully 
occupied  matrix  makes  it  infeasible  in  practice  (Nice  example  for  the  "balance 
of  difficulties",  Sunkel  1984a). 

The  logical  solution,  therefore,  is  to  combine  the  advantages  and  avoid  the 
disadvantages  of  each  method.  See  Lachapelle  (1975,  1977)  for  the  computation 
of  deflections  of  the  vertical,  height  anomalies,  and  gravity  anomalies.  Sunkel 
formulates  the  "ultimate"  combination  as  the  combination  of  LSC  with  integral 
formulae,  high  degree  earth  model,  and  the  best  possible  topographic 
information. 

Turning  to  the  currently  used  procedures  for  the  determination  of  the 
disturbing  potential  above  the  surface  of  the  earth,  the  problem  is  usually 
treated  in  two  steps:  a)  prediction  of  gravity  anomalies  at  the  surface  or  at 
the  sea  (zero)  level,  in  most  cases,  in  the  form  of  mean  anomalies  of  various 
block  sizes;  b)  application  of  one  of  the  boundary  value  solutions  for  the 
computation  of  the  disturbing  potential.  Gravity  anomalies  are  used,  not  only 
because  they  are  the  most  frequently  available  form  of  gravity  information, 
but  also  for  the  reason  that  the  boundary  value  solutions  are  based  on 
gravity  anomalies.  Concerning  the  prediction  of  gravity  anomalies  at  the 
earth’s  surface,  we  are  reminded  of  two  facts  from  the  literature  and 
experience.  One  is  that  the  quality  of  prediction  (interpolation)  depends  on 
the  variance  and  the  correlation  length  of  the  free-air  anomaly  covariance 
function,  in  addition  to  the  data  density.  The  second  fact  is  that  free-air 
anomalies  are  strongly  correlated  with  elevation.  In  areas  of  rugged 
topography,  large  variance  and  small  correlation  length  can  be  expected.  This 
situation  will  result  in  unreliable  interpolated  values.  The  data  should  be 
reduced  to  decrease  the  variance  and  increase'  the  correlation  length.  In 
addition,  the  gravity  data  in  areas  of  rough  topography  are  usually  scarce; 
this  further  aggravates  the  situation.  The  procedure  in  mountainous  areas, 
therefore,  will  be:  removal  of  the  effect  of  topography,  (reduction) 

interpolation  along  the  smoothed  surface  of  the  reduced  data  using  LSC  and 
restoration  of  the  removed  topographical  effect.  Two  of  several  prediction 
techniques  are  recommended  for  the  mean  anomalies  at  the  earth's  surface. 


The  first  is  least  square  collocation  with  the  parameter  model.  The  iterative 
procedure  is  based  on  terrain  corrected  (TC)  free-air  anomalies;  details  are  in 
Svinkel  and  Kraiger  (1983).  It  was  tested  in  Austria  both  in  flat  and  rugged 
areas,  and  it  is  considered  the  "most  practicable"  for  free-air  anomaly 
prediction  along  the  earth’s  surface.  A  fast  algorithm  by  Sideris  (1984) 
computes  terrain  corrections  on  a  grid  using  gridded  topographical  data  and 
FPT  technique.  The  second  technique  is  LSC  and  topographic-isostatic 
reduction  (TIR).  This  approach  is  well  known  and  is  frequently  used.  This 
method  also  assumes  a  compensation  model  and  an  average  density  value, 
which  are  wrong  in  most  cases.  The  differences  between  the  two  methods  are 
as  follows:  In  the  first  approach  the  effect  of  the  topography  is  removed 
using  a  best  possible  estimate  for  the  average  density  in  terms  of  the 
parameter  b  (Bouger  factor),  and  in  the  second  approach  a  standard  density 
is  used;  in  method  one  the  trend  is  optimally  estimated,  and  in  method  two  the 
isostatic  compensation  (a  user  selected  model)  takes  over  the  parameter  model 
of  the  first  method.  It  is  well  known  that  the  actual  density  compensation  is 
largely  different  from  an  ideal  model,  especially  in  small  local  areas.  Lastly, 
the  first  technique  can  process  only  gravity  anomalies  and  requires  an 
iteration  process,  and  the  second  technique  can  process  any  type  of  data  with 
no  iteration  required.  Both  approaches  have  advantages  and  drawbacks;  the 
first  method  is  considered  equal  or  superior  for  most  circumstances,  except 
when  combination  of  heterogeneous  data  is  required. 

The  Stokes-Pizetti  procedure  for  the  determination  of  gravity  field 
quantities  requires  gravity  anomalies  at  zero  level.  The  analytical  continuation 
of  all  types  of  anomalies  to  the  zero  or  sea  level  can  be  performed  by  various 
methods.  Several  of  these  are  discussed  in  the  report.  In  the  opinion  of  the 
author  of  the  report  (Sunkel  1984),  the  LSC  is  the  ideal  tool  for  analytical 
continuation  of  anomalies,  or  any  other  gravity  field  quantity,  to  zero  level. 
Used  in  combination  with  the  quasi-isostatic  anomaly  of  previously  described 
approach  one,  or  applied  to  the  iBostatic  anomaly  of  approach  two,  the  results 
will  be  analytically  continued  quasi-isostatic  or  analytically  continued  isostatic 
anomalies  at  zero  level,  respectively.  Both  of  these  can  be  directly  used  in 
the  Stokes-Pizetti  integral  to  calculate  the  disturbing  potential  outside  the 
earth’s  surface  (after  adding  the  indirect  effect).  Any  other  gravity  field 
quantity  can  be  obtained  from  the  disturbing  potential  by  applying  the 
corresponding  linear  functional  to  the  disturbing  potential  formula. 


In  the  chapter  "Error  Considerations",  the  errors  entering  into  the 
prediction  of  the  disturbance  vector  are  discussed  in  detail.  The  errors  are 
divided  into  four  major  groups,  each  containing  the  relevant  error  sources. 

1. )  Data  reduction  errors 

a)  errors  due  to  terrain  correction  compensation 

b)  errors  due  to  the  DTM  sampling  rate 

c)  errors  due  to  errors  in  the  DTM  data 

d)  errors  due  to  the  computation  of  TIRs 

e)  errors  due  to  approximation  of  the  analytical  continuation 

2. )  Mean  anomaly  errors 

a)  errors  due  to  data  density  and  distribution 

b)  errors  due  to  data  errors 

3. )  Representation  and  errors  due  to  mean  anomaly  errors 

a)  errors  due  to  use  of  mean  anomalies 

b)  errors  due  to  mean  anomaly  errors 

4. )  Errors  of  data  reduction  effects  on  the  estimated  quantities,  breakdown  as 

in  1.) 

All  of  the  above  error  sources  are  discussed  in  detail  and  their  contributions 
to  the  error  budget  of  the  disturbance  vector  analyzed.  Methods  of 
computations  and  numerical  estimates  are  give  for  most  contributing  error 
components. 

Another  study  by  Sunkel  (1981)  estimates  the  accuracies  of  the  gravity 
disturbance  vector  components  at  high  altitudes  from  a  given  set  of  free-air 
mean  anomalies  and  their  rms  errors  at  sea  level. 

The  arrangement  of  the  mean  anomalies  is  symmetrical  with  respect  to  the 
computation  point.  The  innermost  rectangular  zone,  consisting  of  5’  *  5’ 
blocks  of  mean  anomalies,  covers  an  area  of  2*  *  2*,  3*  *  4*,  and  7*  *  9* 
according  to  varying  cases  of  data  distribution.  The  next  zone  of  15*  *  15’ 
mean  anomalies  covers  an  area  varying  in  size  from  6*  *  8*  to  16*  x  18*.  The 
following  rectangular  area  consisting  of  1*  *  1*  mean  values  covers  and  area 
changing  from  26*  *  30*  to  38*  *  42*.  Another  zone  of  5*  x  5*  mean  anomalies 
has  two  variations  in  coverage  area:  50*  *  70*  and  60*  *  100*.  The 
remaining  surface  of  the  sphere  (180*  x  360*)  is  covered  by  5*  x  5*  blocks, 
having  larger  rms  errors.  In  total,  there  are  nine  variations  of  data 
distribution.  The  assumed  rms  errors  for  the  5’  x  5’  blocks  are  *8  or  *10 


mgal;  for  the  15*  «  15'  blocks  *7  and  *8  mgal;  for  the  1*  *  1*  blocks  ‘4  and 
*5  mgal;  and  for  the  5*  *  5*  outer  zones  *3  and  *5  mgal.  In  total,  18  cases 
(according  to  distribution  and  rms  errors)  were  considered  in  the  numerical 
computations. 

Five  altitudes  for  the  computation  point  were  used:  30,000,  40,000,  70,000, 
100,000,  and  200,000  feet. 

Two  methods  were  used  for  the  computation  of  the  disturbance  vector 
components.  These  are  the  integral  formula  for  the  gradient  of  the  disturbing 
potential,  in  terms  of  free-air  anomalies  and  least  squares  collocation.  For  the 
integral  solution,  the  representation  error  was  calculated  in  the  frequency 
domain.  The  results  of  the  two  methods  differ  by  less  than  10%.  It  was 
expected  that  the  results  will  differ  only  slightly  if  the  gravity  coverage  is 
reasonably  good.  This  was  confirmed  by  the  numerical  computations.  An 
optimal  algorithm  was  developed  for  the  collocation  solution  to  take  advantage 
of  the  regularity  and  symmetry  of  the  data  distribution.  From  the  numerical 
investigations  the  following  conclusions  were  made. 

For  radial  components: 

a)  The  error  decreases  rapidly  if  the  number  of  5'  »  5’  mean  anomalies  around 
the  computation  point  increases; 

b)  The  larger  the  5’  *  5’  data  set,  the  smaller  the  representation  error; 

c)  The  effect  of  data  errors  decreases  with  increasing  altitude; 

d)  The  size  of  data  area  of  strong  contribution  increases  with  the  increase  of 
the  prediction  height,  i.e.  if  the  number  of  data  blocks  are  constant,  the 
prediction  error  will  increase  with  altitude. 

For  horizontal  components: 

a)  For  the  same  data,  the  error  is  more  than  twice  as  large  than  the  error  of 
the  radial  component  and  decreases  only  slowly,  if  the  number  of  5'  *  5’ 
anomalies  increases; 

b)  The  larger  the  data  set,  the  smaller  the  representation  error; 

c)  The  influence  of  data  errors  decreases  with  increasing  altitude; 

d)  The  remote  zone  has  a  very  significant  effect  on  the  results. 

The  results  of  the  error  estimates  of  gravity  disturbance  vector  components  at 
the  considered  altitudes  for  the  18  cases  are  tabulated.  The  main  numerical 
results  are:  The  radial  component  can  be  estimated  with  *1  mgal  accuracy  at 


50,000  ft  altitude  (on  the  basis  of  the  stated  data  sets);  to  obtain  the  same 
accuracy  at  30,000  ft,  the  data  accuracy  (especially  the  5’  *  5’  set)  has  to  be 
increased  by  60%.  Regarding  the  horizontal  component,  with  the  best  data 
distribution  *2.3  mgal  accuracy  can  be  achieved  at  30,000  ft  altitude  (0.5"  in 
the  direction  of  the  gravity  vector).  For  the  achievement  of  *1  mgal 
accuracy,  the  block  sizes  must  be  reduced  by  a  factor  of  2  out  to  a  spherical 
distance  of  about  30*,  and  the  overall  data  errors  must  be  reduced  by  30%. 

5. 2. 5. 2.  Upward  Continuation  of  Gravity  Anomalies.  The  upward  continuation 
of  gravity  anomalies  was  recently  investigated  by  Cruz  and  Laskowaky  (1984). 
Three  procedures  were  compared  in  a  numerical  study  with  real  data.  The 
original  gravity  data  consisted  of  about  18,386  free  air  anomaly  values 
irregularly  distributed  over  a  7*  *  9*  area  in  New  Mexico.  The  simple  Bouger 
anomaly  and  the  terrain  correction  were  also  furnished.  Elevation  data  in  the 
form  of  30”  *  30"  grid  point  values  were  also  available.  5'  *  5’  mean 
anomalies  were  developed  in  the  test  area  in  steps  including:  data  thinning, 
terrain  reduction,  tailoring  of  covariance  function,  and  usage  of  only  ten  of 
the  closest  data  points  to  the  point  of  computation  in  the  collocation 
prediction.  "Collocation  from  the  ten  closest  points"  was  used  to  reduce  the 
size  of  the  computation  load.  Mean  elevations  for  5'  *  5'  blocks  were  also 
computed  in  the  7*  *  9*  test  area  from  the  30"  *  30"  grid  data. 

The  given  mean  anomalies  are  referred  to  the  earth's  surface,  the  upward 
continuation  of  these  anomalies  without  reductions  would  be  possible  by 
collocation,  or  by  a  Bjerhammer-type  solution;  however,  these  were  not  feasible 
due  to  large  matrix  inversion  requirements.  If  the  discrete  estimation 
procedures  are  ruled  out,  the  use  of  upward  continuation  of  continuous 
anomaly  functions  require  that  the  data  be  continued  downward  to  the  level 
surface  (sphere).  To  avoid  the  "major  difficulties"  associated  with  this 
reduction,  the  investigators  designed  and  tested  a  procedure  named  "the 
indirect  method".  In  this  method  the  gravity  anomaly  field  is  separated  into 
three  frequency  ranges.  The  high  frequency  part,  the  effect  of  the  shallow 
topography,  is  modeled  by  directly  integrating  the  gravitational  effects 
without  any  reduction  to  a  level  surface  (by  prism  integration).  From  the 
field  left  after  the  removal  of  the  effect  of  the  topography,  the  low  frequency 
component  was  also  removed  by  using  the  Rapp  180  *  180  (1981)  field.  This 
eliminated  truncation  errors  due  to  the  neglect  of  remote  zone  data.  In  the 
removal  process  it  was  assumed  that  the  data  points  are  on  a  common  level 
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surface  instead  of  their  actual  vertical  position.  The  justification  for  this 
assumption  is  that  the  vertical  gradient  of  a  180-field  is  "expected  to  be 
small"  and  that  it  was  necessary  to  keep  the  computation  time  reasonable. 
After  the  removal  of  the  high  and  low  frequency  fields,  the  residual  medium 
field,  much  smoother  than  the  original,  however,  it  is  still  located  on  the 
earth’s  surface.  To  improve  the  situation,  an  approximative  "reduction"  is 
used.  This  is  the  terrain  correction  of  an  expansion  of  the  topography  to 
degree  180,  implicitly  applied  to  the  residual  anomalies.  Detailed  discussion  of 
this  procedure  is  in  Moritz  (1966).  After  all  these  operations,  the  final 
position  of  the  level  surface,  to  which  the  data  are  "reduced",  is  unknown  and 

assumed  to  coincide  with  the  mean  elevation  of  the  topography  in  the  area  of 

upward  continuation.  This  affects  the  upward  continuation  distance  and  the 
results.  The  sensitivity  of  profile  anomalies  to  1.5  km.  change  in  upward 

continuation  distance  (H0)  for  various  heights  are  tabulated  in  the  report. 
The  percentage  ratios  of  rms  change  in  profile  anomalies,  when  K0  is  reduced 
by  1.5  km,  to  the  rms  values  of  profile  anomalies  at  height  Ho  are:  3.8%,  6.8%, 
and  9.1%  at  ho  =  30km,  10km,  and  5km,  respectively. 

The  final  value  of  the  upward  continued  anomaly  is  the  sum  of  the  three 
components  discussed  above.  The  results  of  the  "indirect  method"  was 
compared  with  the  results  obtained  by  the  two  other  procedures.  These  were 
the  "direct  Poisson"  integration  of  the  original  uncorrected  surface  anomalies 
assuming  that  they  are  on  a  level-  surface;  and  the  second,  the  Poisson 

integration  of  terrain  corrected  surface  anomalies  (Faye  anomalies). 

Test  profiles  were  computed  by  the  three  procedures  at  30,  10,  and  5km 
altitudes.  The  results  show  that  the  Poisson  integration  of  terrain 

uncorrected  data  are  too  low  by  0.5,  0.6,  and  0.7  mgal  at  30,  10,  and  5km, 
respectively,  compared  with  the  results  of  the  Poisson  integration  of  terrain 
corrected  data  (this  bias  is  the  upward  continued  terrain  correction).  No  bias 
was  found  between  the  two  terrain  corrected  methods.  Finally,  the  standard 
deviation  of  the  differences  of  all  three  methods  are  the  order  of  (0.5,  0.6, 
1.3)  mgal  at  (30,  10,  5)  km  altitudes. 

Test  computations  were  made  by  the  use  of  FFT  technique  and  compared 
to  the  Poisson  integration  results.  At  30  and  10km  altitudes,  the  agreement 
between  the  two  methods  was  0.1  and  0.3  mgal  respectively. 
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5.2.6.  Airborne  Gradiometry  Improvement  of  the  Local  and  Regional  Gravity 
Field) 

As  it  is  shown  in  Table  2,  section  5.1,  the  high  and  very  high  frequency 
parts  of  the  gravity  field  can  be  obtained  today  almost  exclusively  by  point 
gravity  anomalies.  5*  *  5’  mean  gravity  anomalies  and  deflections  of  the 
vertical  contributed  to  the  high  frequency  field;  however,  the  distribution  of 
these  data  is  mostly  irregular.  The  coverage  of  the  deflections  is  local  and 
the  density  of  data  in  most  cases  is  medium.  Height  data  can  furnish  the 
effect  of  the  visible  topography  and  can  furnish  high  and  very  high  spectral 
resolution  in  combination  with  airborne  gradiometer  data.  An  economically 
feasible  and  efficient  technique  for  observation  of  the  high  frequency  field 
with  high  data  density,  regional  coverage,  regular  data  distribution  and  high 
data  accuracy  will  be  the  airborne  gradiometry.  The  research  and 
development  of  a  feasible  moving  base  gradiometer  system  is  continuing  during 
the  past  15  -  20  years.  The  literature  of  these  developments  is  quite 

extensive,  well-known,  and  therefore  will  not  be  discussed  here.  Nevertheless, 
an  operational  system  is  expected  to  be  completed  soon  and  actual  flight 
testing  can  begin.  Survey  accuracy,  simulation  and  geodetic  application 
Btudies  have  been  made  as  early  as  1975  and  Btill  continuing  today.  Some  of 
these  studies  are:  Moritz  (1975),  Schwarz  (1976,  1977),  White  (1980),  Jekeli 
(1983,  1984),  TASC  (1984). 

A  recent  study  on  achievable  accuracies  for  geodetic  use  by  an  airborne 
gravity  gradiometer  system  is  by  Jekeli  (1983).  The  accuracy  of  the  point 
gravity  vector  on  the  earth’s  surface,  estimated  from  the  five  independent 
components  of  the  second  order  gradient  tensor,  observed  at  altitude,  was 
evaluated  by  a  numerical  analysis.  For  the  method  of  the  analysis  the 
least-squares  collocation  was  chosen  as  the  ideally  suited  optimal  process.  In 
addition  to  the  advantages  of  the  method,  especially  the  immediate  by-product, 
the  prediction  error,  and  the  disadvantages  are  also  discussed.  Among  other 
problem  areas,  Jekeli  points  out  that  the  covariance  model  may  not  be 
adequate  because  the  gravity  field  is  not  known  in  sufficient  detail  for  the 
development  of  a  model  at  very  high  frequency.  It  is  stated,  however,  that 
the  error  analysis  based  on  least-square  collocation,  in  spite  of  some  negative 
aspects,  should  provide  "a  fair  indication"  of  the  performance  of  the 
gradiometer  system.  The  analysis  was  formulated  in  both  the  space  and 
frequency  domains  with  respect  to  either  one  or  two  horizontal  dimensions. 


Approximations  for  a  manageable  analysis  were:  flat  earth  surface,  infinitely 
long  tracks,  continuous  data,  and  infinitely  many  parallel  tracks.  The  analysis 
was  carried  out  with  a  cutoff  wavelength  of  500km  implying  that  the  error 
estimates  refer  to  a  high  frequency  reference  field  corresponding  to  degree 
n  >  80  in  harmonic  expansion.  Other  parametric  values  were:  aircraft  speed, 
300  km/hr;  altitude,  600  m;  gradiometer  white  noise  35Ea/H2,  corresponding  to 
1.9E  standard  error  for  10  seconds  averaging  time;  and  track  spacing,  5km. 

From  the  results  of  the  analyses,  it  is  concluded  that  the  accuracy  of 
gravity  estimated  from  airborne  gradiometer  data  depends  on  the  density  and 
coverage  of  the  survey  tracks.  With  track  spacing  of  >  5  km  and 

unidirectional-track  survey,  an  accuracy  of  1  mgal  cannot  be  achieved.  A 
bidirectional-track  survey  substantially  improves  the  accuracy  of  the  estimated 
gravity.  Track  spacing  and  survey  altitude  are  significant  contributors  to 
total  error.  A  diagram  depicts  accuracies  of  the  gravity  disturbance  (fig)  and 
components  of  the  deflection  of  the  vertical  ((,tj),  estimated  from  gradient 
components  observed  at  altitude  along  bidirectional  tracks,  in  function  of 
track  spacing.  A  second  diagram  shows  the  Bame  in  function  of  survey  height 
above  the  earth’s  surface.  The  accuracy  of  fig  at  5  km  track  spacing  is  about 
0.4  mgal,  at  10  km  spacing  it  is  about  0.7  mgal,  and  at  15  km  track  spacing  it 
is  about  1  mgal.  the  deflection  of  the  vertical  components  accuracies  are  0106, 
Oil,  and  0122  (arc  seconds)  respectively. 

From  the  survey  height  diagram,  at  600  m  the  accuracies  are  as  above. 
At  2000  m,  fig  has  an  error  of  0.6  mgal,  and  at  4000  m  1  mgal.  The  deflection 
components  at  2000  m  have  an  error  of  0109  and  at  4000  m  01 14  (arc 
seconds). 

It  is  noted  that  the  above  figures  will  be  spoiled  somewhat  by  the  used 
approximations,  and  they  are  relative  to  the  reference  field.  If  the  reference 
model  would  be  error-free,  the  above  accuracy  figures  would  represent  the 
absolute  accuracy.  The  current  models  have  quite  large  errors,  especially  at 
the  shorter  wavelengths;  therefore,  the  attainment  of  the  1  mgal  accuracy  goal 
cannot  be  firmly  substantiated. 

Another  analysis  by  Jekeli  (1984)  concerns  the  techniques  for  processing 
of  gradiometry  data.  An  airborne  gravity  gradiometer  system  will  produce  a 


very  large  amount  of  data,  especially  if  the  1  mgal  accuracy  for  an  area 
survey  is  the  desired  goal  and  all  the  six  gradients  at  each  point  are  used. 
In  a  grid  of  300  km  with  5  km  point  density,  there  will  be  21600 


observations. 

If  the  optimal  method,  the  least-squares  collocation,  iB  used,  the  inversion 
of  the  autocovariance  matrix  (21600  *  21600)  generally  would  require  O(1013) 
operations.  Therefore,  an  unmodified  application  of  the  least-squares 
collocation  is  not  feasible,  under  the  present  and  foreseeable  future  computing 
capabilities.  As  in  many  other  applications  of  this  technique,  the 
computational  problem  can  be  improved  without  hurting  the  optimality  of  the 
solution.  An  example  is  applying  certain  restrictions  to  the  pattern  of  the 
observations. 

If  single  observations  are  uniformly  spaced  on  a  straight  line,  the 
auto-covariance  matrix  will  have  a  Toeplitz  structure.  Such  a  simple  Toeplitz 
matrix  can  be  inverted  in  0(n2)  operations,  instead  of  0(n3)  operations 
required  for  the  inversion  of  a  general  auto-covariance  matrix  (n  =  number  of 
observations).  Using  Past  Fourier  Transform,  this  can  be  reduced  to  0[n(log 
n)2]  operations  according  to  the  claim  of  Bitmead  and  Anderson  (1980). 

If  the  measurement  pattern  is  a  grid  of  parallel  tracks  in  both  directions, 
with  uniform  point  spacing,  and  with  one  observation  at  each  point,  the 
auto-covariance  matrix  is  a  "Toeplitz  -  block  Toeplitz  matrix".  If  the  number 
of  observations  at  each  point  is  greater  than  one,  we  have  a  "blocked  Toeplitz 
-  block  Toeplitz  matrix".  An  algorithm  for  the  inversion  of  a  block  Toeplitz 
matrix  in  0(k3n2)  operations  is  given  in  Wiggins  and  Robinson  (1965).  The 
dimension  of  the  block  is  k  *  k.  This  can  be  reduced  to  O[k2n(log  n)2  ] 
(Bitmead  and  Anderson  1980).  An  algorithm  by  Wax  and  Kailath  (1983)  requires 
0(k3n2)  operations  to  invert  a  Toeplitz  block  Toeplitz  matrix.  This  may  be 
transformed  into  one  requiring  0(n3k2)  operations,  and  it  is  beneficial  if 
n  <<  k.  Applying  the  Bitmead  and  Anderson  (1980)  algorithm  to  the  example  of 
a  300  km  grid  with  5  km  uniformly  spaced  points  and  six  gradients  at  each 
point,  the  inversion  will  require  0(  2  *  107)  operations,  which  is  feasible. 

If  rectangular  coordinates  are  replaced  by  cylindrical  polar  coordinates, 
then  the  covariance  matrix  for  uniformly  spaced  single  observations  on  a 
circular  track,  in  addition  to  its  Toeplitz  structure,  is  also  circulant.  Its 
Fourier  transform  is  a  diagonal  matrix;  its  inversion  requires  0(n  log  n) 
operations.  If  we  have  multiple  uniformly  spaced  concentric  tracks  with 
multiple  observations  at  each  point,  then  we  have  the  case  of  blocked 
circulant-block  Toeplitz  matrix.  The  transformation  of  rectangular  pattern  to 
circular  pattern  is  an  additional  processing  of  data. 
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The  second  type  of  technique  discussed  by  Jekeli  (1984)  under  the  title  of 
"Non-optimal  Estimation"  is  reduction  of  the  number  of  observations  by 
averaging.  The  Analytic  Sciences  Corporation  (TASC)  under  contract  to  AFGL 
is  working  on  a  method  where  the  survey  area  around  the  computation  point 
is  divided  into  zones.  The  size  of  the  zones  increases  with  the  distance  from 
the  point.  On  the  basis  of  the  average  observations  in  each  zone, 
least-squares  collocation  is  used  for  the  estimation  of  the  gravity  vector.  It 
is  expected  that  the  differences  between  this  and  the  optimal  solution  will  be 
small.  Simulation  tests  and  comparative  analyses  lie  in  the  future. 

Another  approach  in  the  "non-optimal"  category  is  the  use  of  integral 
formulae  derived  from  the  second  derivatives  of  the  disturbing  potential.  The 
formulae  are  given  in  Jekeli  (1984).  The  formulae  are  exact  except  for  the 
planar  approximation  of  the  earth  (this  is  usually  negligible  considering  the 
extent  of  the  survey).  The  error  sources  are  the  discrete  and  noisy 
observations  and  the  finite  extent  of  the  data  coverage.  The  effect  of  this 
can  be  improved  by  the  modification  of  integration  kernel  (Jekeli,  1982). 
Using  the  integration  formulae  the  disturbance  vector  would  be  at  altitude, 
and  a  simple  least-squares  collocation  procedure  would  compute  the 
corresponding  surface  values. 

In  estimating  anomalous  gravity  or  gradients,  two  or  three  digit  accuracy 
in  the  final  values  are  generally  satisfactory.  Therefore,  under  certain 
circumstances,  a  large  covariance  matrix  can  be  inverted  by  an  approximate 
method  which  is  computationally  feasible.  The  solution  is  non-optimal,  but  the 
numerical  results  are  indistinguishable  from  the  results  of  the  optimal  solution. 
Jekeli  terms  this  type  of  estimation  as  "virtually  optimal  estimation".  An 
example  is  TASC’s  GEOFAST  algorithm  (Tait  1979).  A  class  of  methods  for 
approximation  of  the  inverse  of  a  large  matrix  is  to  improve  sequentially  on 
initial  estimates  of  the  inverse  by  iteration.  These  techniques,  called  the  class 
of  relaxation  methods,  converge  satisfactorily;  a  good  initial  estimate  will  save 
computer  time.  The  techniques  are  described  in  Schwarz  et.  al  (1973).  Jekeli 
investigated  the  conjugate  gradient  method  for  the  inversion  of  the  covariance 
matrix.  From  the  several  options  he  chose  a  separable  matrix  approximation 
suggested  by  Tait  (1979)  for  the  initial  estimate. 

A  numerical  analysis  was  performed  for  the  purpose  of  comparison  of 
several  data  processing  techniques.  The  assumed  measurement  grid  was  a 
regular  one  with  data  points  5  km  apart  in  each  direction.  The  extent  of  the 
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grid  was  variable  from  50  km  to  300  km,  and  the  altitude  above  the  ground 
was  600  m.  The  assumed  accuracy  of  the  gradient  observation  was  2E 
(uncorrelated  random  error).  The  average  accuracies  of  estimating  the 
residual  gravity  disturbances,  using  the  (1)  integral  formulas,  (2)  integral 
formulas  with  modified  kernel,  (3)  conjugate  gradient  method  with  iterations, 
and  (4)  the  least  squares  collocation  method,  are  compared  in  one  tabulation 
and  five  diagrams.  The  overall  conclusions  are  that  the  least-squares 
collocation  is  the  optimal  method  in  terms  of  accuracy,  but  requires 
"inordinate"  amount  of  computation  "in  the  most  general  case".  The  integral 
formulae  may  be  sufficient  in  some  applications,  but  not  for  the  achievement  of 
1  mgal  accuracy. 

The  conjugate  gradient  technique  for  the  inversion  of  the  covariance 
matrix  of  the  least-squares  collocation  looks  promising  regarding  reductions  in 
computations.  It  was  experienced  that  convergence  to  1  or  2  digits  is 
achieved  in  much  fewer  iterations  than  the  theoretically  required  n  iterations, 
where  n  is  the  size  of  the  matrix. 

An  earlier  study  regarding  achievable  accuracies  from  airborne 
gradiometry  was  performed  by  Schwarz,  K.P.  (1976),  and  a  simulation  study 
also  by  Schwarz  (1977). 

The  accuracy  study  investigates  the  achievable  accuracies  of  mean  gravity 
anomalies  for  various  block  sizes  from  5’  *  5’  to  5*  *  5*  from  observed  first 
and  second  order  gradients  at  flying  altitudes.  Gravity  anomalies  are  derived 
by  integrating  the  first  order  gradients  observed  by  accelerometers 
simultaneously  with  the  gradiomoter  system.  In  other  words,  this  concept  can 
be  interpreted  as  a  combined  set  of  measurements  where  the  integrated  first 
order  gradients  provide  the  reference  field  (so  to  speak)  for  the  fine 
structure  observed  by  the  gradiometer  system.  The  theoretical  background  of 
separation  of  gravitation  and  inertia  and  the  derivation  of  related  mathematical 
formulae  are  given  in  several  publications  of  Moritz  (1967,  1971,  1975).  In 
Moritz  (1975)  and  Schwarz  (1976),  a  system  of  linear  second  order  differential 
equations  are  given  interrelating  the  gravity  disturbance  vector  with  the 
observed  second  order  gradients  and  accelerometer  outputs  (gravitational  plus 
inertial  force  combined).  For  the  combination  of  the  measured  second  order 
gradients  and  integrated  values,  their  downward  continuation  from  the  flight 
altitude  to  the  surface  of  the  earth  and  for  the  computation  of  mean  values, 
the  method  of  least  squares  collocation  was  used.  All  of  the  above  steps  were 


performed  by  appropriate  selection  and  changes  of  the  covariance  function. 
In  the  particular  case  of  Schwarz's  study,  the  combination  with  satellite 
altimeter  data  was  also  included  in  the  algorithm. 

Due  to  the  decisive  role  of  the  covariance  function  in  the  collocation 
method,  after  a  careful  study,  two  models:  Hirvonen’s  (1962),  derived  from 
local  data;  and  Kaula's  (1966),  derived  from  local  and  regional  data,  were 
selected.  Using  the  characteristic  parameters  of  both  types  of  functions 
(Hirvonen’s  called  "local"  and  Kaula's  called  "regional"),  the  covariance 
functions  of  the  gravity  anomalies,  second  order  radial  gradients,  and  geoid 
undulations  have  been  computed  for  0  km  (sea  level)  and  10  km  (flight) 
altitudes. 

The  measurements  used  are,  as  previously  mentioned,  gradiometer  data  in 
a  (♦,A,r)  -  system  and  Ag  values  from  integration  along  the  profile.  Mean 
values  at  10  second  intervals  were  used  as  gradiometer  measurements.  For 
the  consideration  of  the  measurement  errors,  the  covariance  matrix  of  the 
measurements  (Cxx)  was  augmented  by  the  covariance  matrix  of  the  noise.  To 
obtain  information  on  the  sensitivity  of  the  estimation  to  the  errors  of  the 
different  measurements,  the  standard  errors  of  each  type  of  measurements 
have  been  varied  separately.  The  first  values  estimated  were  the  Ag  values  in 
any  arbitrary  position  between  two  flight  profiles.  This  can  be  done  by 
combination  of  various  data  of  the  respective  flight  profiles.  The  following 
combinations  were  considered:  (1)  Ag  values  interpolated  between  the  profiles, 
(2)  use  of  the  three  first  order  gradients  condensed  from  the  five  gradiometer 
measurements,  (3)  five  second  order  gradients  alone,  and  (4)  a  combination  of 
Ag  and  different  second  order  gradients.  The  numerical  problem  is  the 
inversion  of  the  covariance  matrix  of  the  measurements.  After  some  test 
computations  it  was  found  that  an  interpolation  point  between  the  flight 
profiles  can  be  satisfactorily  determined  by  using  a  range  of  observation 
points  along  the  respective  flight  profiles.  The  length  of  the  range  is  about 
three  times  the  distance  from  the  estimated  point  to  the  profiles.  It  was 
determined  that  the  standard  error  decreases  only  slightly  if  more  points  are 
used.  The  interpolated  point  should  be  an  equal  distance  between  the  two 
profiles.  Computations  were  performed  using  both  the  "local"  and  "regional" 


covariance  functions.  For  the  purpose  of  the  accuracy  study,  only  the  error 
covariance  matrix  was  computed. 

The  downward  continuation  from  flight  altitude  to  zero  level  was 
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performed  by  the  selection  of  a  covariance  function  that  can  be  continued 
analytically  to  the  sea  level,  so  that  the  downward  continuation  became  a 
problem  of  spatial  interpolation.  Like  in  the  planar  interpolation  at  the  flight 
altitude,  the  accuracy  of  the  interpolated  point  depends  on  the  distance  to  the 
observation  points.  In  both  cases  (planar  and  special  interpolation)  the  best 
data  combination  is  Ag-values  combined  with  gradiometer  data.  The  optimal 
solution  is  the  combination  of  Ag  and  five  independent  gradients.  For 
east-west  profiles,  the  combination  of  Ag,  Tr$,  and  is  very  close  to  the 

best  combination. 

To  asses  the  influence  of  the  measurement  errors,  two  accuracy  diagrams 
are  given  in  the  report.  One  shows  the  standard  error  of  the  estimation  in 
the  function  of  the  standard  error  of  the  observed  gravity  anomalies  and 
errorless  gradient  measurements;  the  second  diagram  shows  the  estimation 
error  in  the  function  of  the  standard  error  of  the  gradient  measurements. 
Both  error  diagrams  give  the  estimation  errors  for  profile  spacings  of  0.2*, 
0.5*,  0.7*,  and  1.0*.  For  the  computation  of  mean  gravity  anomalies  for 
various  block  sizes,  the  basic  least-squares  collocation  formula  used  for  the 
point  anomalies  was  modified.  For  the  derivation  of  the  required 
cross-covariance  function  between  point  and  mean  anomalies  smoothing 
operators  were  computed  for  blocks  of  different  BizeB  according  to  Meissl 
(1971).  Using  different  assumptions  for  the  measurement  errors  and  various 
profile  spacings  mean  gravity  anomaly  accuracies  were  computed  for  several 
block  sizes. 

The  conclusions  of  Schwarz’s  accuracy  study  are  summarized  as  follows: 

1.  The  accuracy  of  an  interpolated  anomaly  depends  on  its  distance  from  the 

nearest  profile.  Therefore,  the  dense  spacing  of  the  observation  profiles 

is  more  desirable  than  cross  profiles.  Cross  profiles  should  be  used  for 
updating  only. 

2.  A  high  density  of  observed  points  along  a  profile  is  necessary  only  for 
the  integration  of  first  order  gradients.  For  interpolation  of  Ag  values,  a 
density  between  1/4  and  1/2  of  the  profile  spacing  is  sufficient. 

3.  Optimal  data  combination  for  profile  spacing  larger  than  0.3*  should 

include:  Ag,  Tr4,  Trx,  T*4,  Txx»  for  spacings  smaller  thar  0.3*,  Ag,  Tr4, 

Trx>  and  Trr  are  sufficient. 

4.  The  effect  of  measurement  errors  is  larger  for  small  profile  spacing. 

5.  The  combination  with  altimeter  data  is  significant  for  Large  (5*  *  5*)  mean 


values.  The  contribution  to  1*  *  1*  blocks  is  marginal. 

6.  The  attainable  accuracies  under  the  assumption  made  about  the  covariance 
functions  were: 

a)  with  1*  profile  spacing,  an  accuracy  of  *3  mgal  standard  error  for 
5*  ■  5*  blocks  and  ‘5  mgal  for  1*  *  1*  blocks  can  be  achieved  (for  the 
5*  >  5*  blocks,  only  with  the  combination  of  altimeter  data). 

b)  with  0.3*  profile  spacing  for  15’  *  15’  and  5’  *  5’  blocks,  *3  mgal  can 
be  obtained. 

In  the  "Simulation  Study  of  Airborne  Gradiometry”  (Schwarz  1977) 
simulated  gravity  anomalies  and  second  order  gradients  are  used  for  the 
recovery  of  the  anomalies  at  ground  level  from  the  simulated  measurements  at 
flight  level.  The  purpose  of  the  study  was  to  obtain  an  algorithm  for 
airborne  gradiometry,  and  to  design  a  method  for  efficient  handling  of  very 
large  amounts  of  data  expected  from  the  measuring  Bystem.  For  covering  an 
area  of  20*  *  25*  with  profiles  spaced  at  1*  and  assuming  250  observations 
per  profile  degree,  there  will  be  130,000  measurements;  if  20’  profile  spacing  is 
used,  the  number  of  observations  will  be  390,000. 

The  study  consists  of  four  distinct  steps:  (1)  gravity  and  gradiometry 
data  are  generated  at  ground  and  flight  level,  (2)  the  flight  level  data  are 
corrupted  by  an  error  model,  (3)  the  ground  level  gravity  anomalies  are 
estimated  from  the  flight  level  data,  and  (4)  estimated  and  generated  ground 
level  data  are  compared. 

The  gravity  field  data  was  generated  from  a  grid  of  mass  points  of  10’ 
spacing  at  40  km  depth  for  an  area  of  10*  ■  15*.  A  positive  or  negative  mass 
of  constant  size,  or  a  zero  mass,  was  selected  according  to  a  normal 
distribution  (the  data  have  been  generated  by  DMAAC,  St.  Louis). 

For  the  estimation  of  gravity  anomalies  at  ground  level  from  the  simulated 
observations  at  flight  level  (10  km),  the  formula  of  least-squares  collocation 
was  used.  According  to  the  spacing  of  the  profiles  and  the  size  of  the  mean 
anomalies  to  be  estimated  an  optimal  point  configuration  was  chosen.  It  was 
assumed  that  the  observation  profiles  are  parallel  and  the  observation  points 
are  at  constant  intervals  (At  500  knots  speed  and  10  seconds  integration  time, 
the  observation  rate  was  a  uniform  2.6  km).  With  the  above  conditions  a 
"moving  operator"  was  computed  (from  the  cross-covariances  between  signal 
and  observation  and  from  the  auto-covariance  matrix  of  the  observations). 
This  unchanged  operator  was  moved  along  the  flight  profiles  from  one  set  of 


data  points  to  the  next,  according  to  the  pattern  described  previously  in 
Schwarz  (1976),  using  all  available  observations. 

The  covariance  function  used  in  this  study  was  that  used  in  Schwarz 
1 1976).  Its  parameters  are  different  Bomewhat  from  the  covariance  function 

directly  derived  from  the  simulated  data.  This  was  also  used  for  the 

determination  of  the  effects  of  assumption  errors  in  the  covariance  on  the 
results  of  the  estimation.  For  the  statistical  error  models  "normal  models” 
(numbers  taken  from  a  normal  distribution),  and  Markov  sequences  of  first 
and  second  order  were  used.  The  difference  between  the  normal  and  Markov 
model  is  that  normal  deviates  are  uncorrelated,  elements  of  Markov  sequences 
are  not.  To  illustrate  the  accuracy  of  the  estimated  point  anomaly,  a  3*  long 
"interpolation  profile"  on  the  ground  level  obtained  from  "flight  profiles" 
spaced  at  20'  is  compared  to  the  exact  profile  of  the  simulated  field.  The 
normal  error  model  was  used,  and  the  standard  errors  of  the  observed  gravity 
anomalies  and  of  the  second  order  gradients  were  assumed  as  *1  mgal  and  *1E 

respectively  (too  optimistic  assumptions).  The  agreement  between  the  two 

plotted  profiles  is  very  good.  The  standard  error  is  3.3  mgal  for  the  point 
estimation.  If  the  interpolated  profile  is  determined  form  the  same  data  as 
before,  but  is  located  directly  below  one  of  the  flight  profiles,  the  standard 
error  drops  to  1.96  mgal.  It  is  concluded  that  a  measuring  error  of  *1  mgal 
at  flight  altitude  (10  km)  is  amplified  to  about  * 2  mgal  by  downwa  d  contin¬ 
uation,  and  to  about  *3  mgal  when  interpolated  between  two  profiles  20’  apart. 

The  results  for  the  mean  anomalies  are  tabulated  and  compared  to  the 
accuracies  obtained  by  the  accuracy  study  (Schwarz  1976).  With  profile 
spacing  of  20’,  the  accuracies  of  the  15’  x  15’,  30’  x  30’,  and  1*  x  1*  blocks 
agree  well  with  the  values  obtained  by  the  accuracy  study.  The  average  of 
the  differences  between  the  two  sets  of  accuracy  figures  is  about  0.4  mgal. 
For  the  1*  «  1*  blocks  obtained  from  1*  profile  spacing,  the  simulation  study 
gave  *4.6  mgal,  versus  the  *5.6  mgal  obtained  by  the  accuracy  study.  The 
normal  error  model  was  used  for  the  above  computations. 

In  the  conclusions  it  is  pointed  out  that  the  two  studies  (accuracy  and 
simulation)  agree  well  regarding  the  obtained  accuracies,  however,  it  is  also 
pointed  out  that  correlated  errors  in  the  measurements  will  strongly  effect  the 
accuracy  of  the  results.  Depending  on  the  size  of  the  correlations  and  on  the 
variances,  the  mean-square  errors  may  double  as  compared  to  the  uncorrelated 


6.  SUMMARY  AND  RECOMMENDATIONS 


The  principal  requirement  for  the  derivation  of  a  model  which 
satisfactorily  approximates  the  earth's  gravity  field  is  a  uniform,  dense,  and 
global  coverage  of  several  types  of  gravity  data.  The  techniques  to  provide 
the  instrumentations  necessary  for  such  coverage  already  exist  or  are  in  their 
final  stages  of  development.  The  achievement  of  ambitious  goals,  like  NASA’s 
Geopotential  Research  Mission  (GRM),  would  require  additional  substantial 
efforts  in  data  collection  programs  and  data  processing  techniques. 

Current  geopotential  models  are  derived,  generally,  from  the  combination 
of  three  types  of  information:  perturbation  of  satellite  orbits,  satellite 

altimeter  data,  and  terrestrial  mean  anomalies. 

Satellite  orbital  data  represent  the  long  wavelength  part  of  the  gravity 
spectrum,  the  short  wavelength  parts  are  smoothed  out  by  the  attenuation  of 
the  gravity  field  with  the  altitude.  Models  containing  well  distributed  and 
precise  orbital  perturbation  data  determine  the  low  order  coefficients  very 
well  (i.e.  degree  and  order  4)  and  they  are  very  useful  for  satellite  orbit 
computations.  Examples  are  GEM9,  GEM-L2,  GRIM-3B,  and  GRIM-MPI. 

From  satellite  altimetry,  the  geoid  over  the  oceans  is  obtained.  These  data 
are  regularly  distributed  and  could  yield  about  50  km  half  wavelength 
resolution.  If  so,  this  would  give  an  expansion  of  360  degree  and  order, 
provided  the  coverage  extends  over  the  globe.  currently,  1*  *  1*  mean 

anomalies  are  developed  from  the  altimeter  data  and  merged  with  gravity 
anomalies  of  1*  *  1*  obtained  over  land.  One  of  the  two  problems  with  this 
type  of  data,  the  effect  of  averaging,  can  be  reversed  by  the  use  of 
"desmoothing"  factors  (Colombo  1981a).  The  other  problem,  that  of  aliasing, 
remains  after  desmoothing  or  filtering.  In  the  case  of  1*  *  1*  mean  anomalies, 
the  aliasing  errors  are  about  50%  of  the  coefficient  at  the  Nyquist  frequency, 
and  larger  above.  Therefore,  there  is  a  need  to  replace  the  current  mean 
values  with  smoothed  values  free  of  aliasing.  It  is  recommended  to  study  this 
problem  and  to  test  K.P.  Schwarz’s  (1984)  suggestion  "to  use  a  bandpass  filter 
on  the  power  spectrum,  obtain  the  auto-covariance  function  by  Fourier 
transforming  the  results  to  the  space  domain  and  estimate  smoothed  values  in 
the  center  of  the  block  using  thiB  function". 

In  the  combination  of  various  data  sets  for  global  models,  or  using  global 
models  as  reference  fields  for  regional  or  local  gravity  field  approximations, 
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the  proper  weighting  of  the  various  data  set8|  or  models,  would  require  their 
actual  accuracies.  Wenzel  (1982)  proposed  the  use  of  the  error  covariance 
functions  of  various  data  types  for  the  derivation  of  spectral  weighting 
functions  for  use  in  data  combinations.  This  approach  can  be  used  in 
connection  with  any  combination  technique  and  it  is  "most  promising"  (Schwarz 
1984).  It  is  recommended  to  investigate  the  possibility  of  whether  the  spectral 
weighting  functions  of  a  particular  data  set  could  establish  a  link  between  the 
frequency  range  of  the  data  sets  and  the  accuracy  of  the  combined  product, 
or  harmonic  coefficients  of  a  certain  degree  and  order  derived  from  the  data. 

From  the  inspection  of  the  percentage  errors  of  the  coefficients  of  various 
solutions  (e.g.  Rapp  1981),  it  can  be  seen  that  the  coefficients  above  degree  12 
have  larger  errors.  It  is  also  apparent  that  up  to  degree  15,  the  coefficients 
are  influenced  mainly  by  satellite  perturbation  data,  and  above  degree  15  the 
mean  anomalies,  including  altimeter  data,  dominate  the  accuracy  of  the 
combination. 

The  accuracy  of  current  models  are  illustrated  by  the  intercomparison  of 
the  result  obtained  by  various  groups.  It  is  true  that  the  various  groups 
frequently  used  the  same  observed  material,  but  their  combination  and 
adjustment  methods  are  different.  The  comparisons  indicate  only  the  order  of 
magnitude  of  the  errors  in  the  various  solutions.  Generally  the  data  errors, 
lack  of  data  continuity  and  aliasing  are  the  principal  sources  of  error.  For 
details  see  section  2. 


6.1.  Satellite  Programs  for  the  Improvement  of  the  Gravity  Field 

6.1.1.  Satellite  Altimetry 

The  most  productive  program  contributing  to  the  improvement  of  global 
gravity  coverage  was  satellite  altimetry.  GEOS-3  and  SEASAT  results,  their 
contributions  to  the  geodetic  aspects  of  the  program,  and  future  plans  for 
satellite  altimetry  are  discussed  in  detail  in  Section  3.1  of  this  report. 

As  it  was  mentioned  earlier,  the  disadvantage  of  satellite  altimetry  from 
the  point  of  view  of  geopotential  modeling,  is  that  the  data  coverage,  due  to 
its  very  nature,  is  restricted  only  to  the  ocean  areas.  This  requires  the 
combination  with  another  type  of  data  (terrestrial  mean  anomalies).  Therefore, 
the  global  data  homogeneity  is  lost.  Nevertheless,  a  homogeneous,  regular, 
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and  fairly  accurate  coverage  was  possible  over  the  areas  by  this  technique. 

The  results  of  satellite  altimetry  can  be  improved  upon  by  improving  three 
related  factors:  a)  sea  topography  and  tidal  corrections,  b)  accuracy  of  the 
radial  component  of  the  satellite  orbit,  and  c)  precision  of  the  altimeter 
measurements. 

Currently,  two  new  altimeter  missions  are  in  Borne  stages  of  their 
development:  The  GEOSAT  and  the  TOPEX.  The  GEOSAT  of  the  Navy  was 

scheduled  for  the  mid  1980’s  and  the  TOPEX  of  NASA  is  contemplated  for  the 
1990’s.  GEOSAT  will  have  a  SEASAT  type  of  sensor  and  about  18  months 
nominal  mission  time.  The  plans  for  TOPEX  call  for  an  altimeter  accuracy  of  *2 
cm  and  for  a  sea  surface  height  accuracy  along  a  grid  of  *14  cm. 

For  the  improvement  of  the  radial  component  of  the  orbit,  an  improved 
TRANET  Doppler  network  is  planned,  and  the  Jet  Propulsion  Laboratory  (JPL) 
intends  to  develop  a  new  tracking  system  (Series  X)  using  the  GPS  system  for 
the  determination  of  the  tracking  sites'  coordinates. 

6.1.2.  Satellite  to  Satellite  Tracking  and  Satellite  Gradiometry. 

These  techniques,  discussed  in  detail  in  sections  3.2  and  3.3,  are  similar  in 
the  respect  that  both  have  a  capability  to  provide  a  high  density,  regular  and 
global  coverage  of  gravity  data.  The  "global  coverage"  in  both  cases  are 
real,  i.e.  covering  land  and  sea,  pole  to  pole,  without  any  gaps.  The 
frequency  range  for  both  will  be  in  the  medium  field,  corresponding  to  the 
frequencies  of  1*  *  1*  mean  anomalies,  provided  by  satellite  altimetry,  and  by 
1*  *  l*  land  values.  The  possibility,  depending  on  the  survey  configurations, 
exists  for  extension  into  higher  frequency  ranges,  particularly  for  the 
gradiometer  survey.  The  data  accuracies  will  be  superior  to  the  accuracy  of 
the  current  altimeter  derived  data. 

The  results  of  several  accuracy  and  simulation  studies  are  summarized  in 
section  3.2  and  3.3.  The  results  of  some  of  these  are  given  as  follows: 

For  a  low-low  satellite  to  satellite  tracking  mission  of  six  months,  at  an 
altitude  of  160  km  and  a  range-rate  noise  of  *1  micrometer/sec  with  4  sec 
integration  time,  the  accuracy  of  the  derived  mean  anomalies  of  1*  *  1*  is 
estimated  as  *2.3  mgal,  and  of  the  1*  *  1*  geoid  undulations  as  4.3  cm  (Jekeli 
and  Rapp  1980) 

Colombo  (1981b)  from  the  error  analysis  of  the  global  geopotential, 
obtained  from  a  low-low  SST  mission,  under  certain  conditions,  the  following 


accuracies: 


(1)  The  relative  error  in  the  potential  coefficients  could  be  better  than:  1% 
up  to  degree  n  =  130;  better  than  10%  up  to  n  =  210;  and  better  than  50% 
up  to  degree  n  =  270. 

(2)  The  geoid  point  value  accuracy  could  be  better  than  5  cm  rms  in  the  band 
from  3000  km  to  40,000  km  and  better  than  10  cm  in  the  band  from  140  km 
to  3000  km  (both  total  errors). 

For  satellite  gradiometry  two  instruments  are  under  development:  the 
cryogenic  gravity  gradiometer  by  Paik  (1981)  and  "Project  Gradio"  by  Balmino 
et  al  (1984).  The  instrument  precision  goals  are  10-4E  (IE  =  10-9s-2).  An 
operational  system  is  expected  for  the  1990’s. 

Simulation  studies  for  five  or  one  (radial)  component  systems  are  available. 
Some  are  summarized  in  section  3.3.  The  most  recent  is  by  Rapp  (1985).  The 
SST  and  satellite  gradiometry  is  intercompared  with  two  recent  geopotential 
models.  (Rapp  81  and  GEM-L2).  The  SST  mission  improves  current  models  by 
a  factor  of  10.  The  gradiometer  mission  improves  the  SST  results  by  about 
60%  with  the  exception  of  wavelengths  shorter  than  20  km. 

6.2.  Data  Processing  and  Adjustment  Techniques 

The  processing  at  a  large  amount  of  observations  and  the  combination  of 
several  types  of  gravity  data  for  estimation  of  gravity  field  models  is  today 
not  on  the  optimal  level.  It  is  widely  recognized  that  the  least  squares 
collocation  is  the  optimal  method  for  the  solutions  of  many  problems  involving 
functionals  of  the  gravity  potential.  It  is  also  widely  known  that  in  many 
cases  the  advantages  of  this  method  is  balanced  out  by  the  problem  of 
inverting  large  matrices,  the  sizes  of  which,  normally,  equals  the  number  of 
observations..  The  problem  of  the  number  of  observations  will  be  more 
difficult  if  the  data  from  satellite  to  satellite  tracking  and  satellite  gradiometry 
will  be  available. 

The  problem  can  be  improved  with  utilization  of  the  symmetries  of  data 
and  with  some  restrictions  to  the  pattern  of  the  observations.  Colombo’s 
(1981)  algorithms  for  harmonic  analysis  is  an  example  where  the  symmetries  of 
data  and  relations  between  spherical  harmonics  and  Fourier  series  was  used 
for  the  evaluation  of  harmonic  coefficients  (Section  3.5). 
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Hajela  (1984)  implemented  Colombo’s  (1981)  algorithms,  originally  designed 
for  5*  «  5*  anomalies,  for  global  1*  *  1*  anomalies  with  some  necessary 
modifications.  The  optimal  estimation  of  the  coefficients  was  carried  to  degree 
and  order  250.  The  use  of  optimal  estimation  resulted  in  improvement  versus 
Rapp’s  coefficients  obtained  by  a  good  approximative  process  (Section  3.5). 

Jekeli  (1984)  discussed  how  the  computational  load  can  be  reduced  without 
hurting  the  optimal  solution  of  the  least  squares  collocation.  According  to  the 
various  patterns  of  observations,  the  auto-covariance  function  will  have 
different  Toeplitz  structures.  The  number  of  the  required  operations  for  the 
different  Toeplitz  matrices  are  given  in  the  paper  (see  Section  5.2.6.). 

S.C.  Bose  et  al.  (1983),  like  Colombo's  approach,  recommend  the  more 
extensive  exploitation  of  the  grid  structure  of  the  data  and  of  the  structure 
of  the  covariance  matrix  for  the  reduction  of  the  computations.  They  explore 
the  equatoi  .al  and  rotational  symmetries  to  avoid  the  ill-conditioning  by 
crowded  data  in  polar  regions.  It  is  shown  that  at  high  latitudes  thinned-out 
data  samples  can  be  used  (Section  3.5) 

Rapp  (1984)  reviewed  some  of  the  adjustment  techniques  for  the 
combination  of  satellite  and  terrestrial  data  for  the  development  of  harmonic 
coefficients.  For  the  method  based  on  orthogonality  relation  between  anomalies 
and  harmonic  coefficients,  correction  terms  were  derived  for  the  effects  of  the 
assumptions  of  spherical  approximation  and  of  the  zero  elevation  of  the 
anomalies.  The  effects  of  the  elevation  is  small  (1%  -  3%);  the  spherical 

approximation  error  is  large,  10%  at  degree  75  and  31%  at  degree  180.  Some 
alternate  adjustment  techniques  are  discussed  using  the  orthogonality,  with 
anomalies  on  a  bounding  sphere  enclosing  all  the  topography.  It  is 

recommended  to  use  0.5*  *  0.5*  mean  anomalies  globally  when  this  data 

becomes  available  over  significantly  large  areas  (Section  3.5). 

6.3.  Estimation  of  Gravimetric  Quantities  in  Local  and  Regional  Areas. 

The  estimation  of  the  local  or  regional  gravity  field  usually  consist  of  the 
determination  of  some  functionals  of  the  anomalous  gravity  potential  in  a 

particular  area  or  at  selected  points. 

The  quality  of  the  estimated  gravimetric  quantities  will  be  determined  by: 

(1)  the  density  of  the  available  data 


u  .■  w,  ^  *r, ^  w  i  [r>  r»  v  *  ■:  •  ">  - 


(2)  the  extent  of  the  data  coverage 

(3)  the  sensitivity  of  the  function  to  be  estimated  to  the  frequency  of 
the  given  data  set 

(4)  measurement  accuracy. 

The  spectral  sensitivities  of  geoid  height  (N),  gravity  anomaly  (A g)  and  of 
the  second  order  radial  gradient  (Trr),  as  computed  by  Schwarz  (1984)  are 
listed  in  Table  5.1.  The  current  and  future  data  types  are  characterized 
according  to  spectral  resolution,  density,  coverage,  distribution,  and  accuracy 
in  Table  5.2  (after  Schwarz  1984). 

If  a  truncated  series  of  harmonic  coefficients  is  used  as  a  reference  field 

for  some  local  estimation,  for  the  proper  weighting  of  the  two  data  sets,  it  is 

desirable  to  know:  a)  the  part  of  the  total  spectral  power  contained  in  a 

solution  of  degree  and  order  N,  b)  which  error  spectrum  is  associated  with 

the  global  solution?  This  subject  area  should  be  investigated  in  detail 

(Section  5.1). 

The  local  or  regional  gravimetric  quantities  can  be  computed  by  integral 
formulae  or  by  spherical  harmonic  series.  The  two  are  equivalent 

theoretically.  Practically,  however,  they  are  different  due  to  the  differences 
in  input  data  characteristics.  Since  the  theoretical  data  requirements  for  the 
integration  formulae  cannot  be  fully  satisfied,  these  formulae  have  been 
modified  and  used  in  combination  with  other  data  contributing  the  information 
outside  of  the  zone  of  integration.  Spherical  harmonic  expansions  of  the 
potential  to  degree  and  order  180  (Rapp  1978,  1981;  Lerch  et  al  1981), 

expansions  of  the  topography,  of  the  rock  equivalent  topography  and  of  the 
isostatically  compensated  topography  (Rapp  1982)  are  available.  Gridded 
digital  terrain  models,  (DTMs)  are  also  available  for  many  areas  for  smoothing 
the  local  gravity  field,  necessary  especially  in  topographically  rough  areas. 

A  number  of  studies  for  the  estimation  of  local  and  regional  gravity  fields 
are  listed  and  some  are  reviewed  in  detail  in  Section  5.2.1.  The  modifications 
of  the  Stokes  and  Vening-Meinesz  formulae,  their  combination  with  harmonic 
coefficients  and  topographic  data  are  discussed  in  Lachapelle  (1984), 
Lachapelle  (1978),  and  Jekeli  (1981). 

Gravimetric  and  satellite  derived  undulations  are  compared  by  Rapp  and 
Wichiencharoen  (1984)  using  20  Doppler  sites  in  the  U.S.,  ten  in  the  mountains 

and  ten  in  flat  areas.  The  results  of  the  study  are  given  in  Section  5.2.2. 

The  accuracies  of  height  anomalies  and  deflections  of  the  vertical  obtained 
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from  the  combination  of  harmonic  coefficients  and  mean  gravity  anomalies  was 
studied  by  Heck  (1983).  The  integration,  replaced  by  summation,  is  limited  to 
a  cap  around  the  computation  point,  the  outer  zones  are  replaced,  or 
represented,  by  a  geopotential  model.  The  errors  of  various  sources  are 
analyzed  and  it  is  revealed  that  at  the  zeros  of  the  kernel  functions  within 
the  spherical  integrals,  the  error  functions  show  local  minima.  Therefore,  the 
integration  radius  should  be  extended  exactly  to  the  first  zero  of  the  modified 
Stokes  and  Vening  Meinesz  kernel  functions  (Section  5.2.3). 

Various  methods  of  the  consideration  of  terrain  effects  in  the  estimation  of 
local  gravity  field  quantities  and  the  use  of  the  stepwise  collocation  for  the 
estimation  are  discussed  in  two  studies.  The  first,  by  Forsberg  and 
Tscherning  (1981),  is  summarized  in  Section  5.2.4.  A  number  of  observed 
gravity  anomalies  and  deflections  of  the  vertical  were  predicted  in  the 
mountainous  area  of  New  Mexico.  The  original  unchanged  gravity  data  and 
their  corrected  versions  by  several  reduction  methods  were  used  to  illustrate 
the  effect  of  the  reductions.  The  conclusion  of  the  study  is  that  it  is 
possible  to  predict  deflections  of  the  vertical  and  gravity  anomalies  in  areas 
of  rugged  topography  with  accuracies  of  1"  and  3-4  mgal  respectively  from 
anomaly  data  spaced  at  6'  apart,  if  terrain  corrections  are  computed  from 
0.5  *  0.5  minutes  point  elevation  data.  The  other  study,  "The  geoid  of 
Austria"  by  Sunkel  (1983),  determined  the  geoid  from  521  observed  deflections 
of  the  vertical  and  a  20"  x  20"  digital  terrain  model  as  local  data.  In 
addition,  the  Rapp  1981  geopotential  model  and  a  1*  x  1*  mean  digital  terrain 
model  (DTM)  were  used  as  global  information.  The  results  from  a  stepwise 
collocation,  topographic-isostatic  reduction,  and  of  the  geopotential  model  were 
combined  to  obtain  the  final  values.  The  differences  between  height  anomalies 
and  geoid  heights  are  in  very  good  agreement  with  the  results  obtained  from 
Bouger  anomalies  and  topographical  elevation  data  (Section  5. 2. 4. 2). 

Some  other  collocation  studies  for  prediction  of  free-air  anomalies  (Sunkel 
and  Kraiger  1983),  combination  of  collocation  and  geophysical  inversion,  and 
studies  of  terrain  reduction  accuracies  in  the  report  of  Forsberg  (1984)  are 
reviewed  in  Section  5. 2. 4. 3. 


6.4.  Gravity  at  Altitudes 


Several  of  the  more  recent  studies  on  this  subject  are  described  in 
sections  5.2.5. 1  and  5. 2. 5. 2. 

In  Sunkel  (1984a),  four  different  approaches  are  discussed  for  the 
prediction  of  gravity  disturbance  at  altitudes,  with  the  advantages  and 
disadvantages  of  each.  The  recommended  most  feasible  method  is  a 
combination  of  least  squares  collocation,  integral  formulae,  high  degree  earth 
model,  and  the  best  topographic  information  (the  last  item  is  crucial  in  areas 
of  rugged  topography).  A  detailed  error  analysis  is  given  for  each  step  of 
the  process,  with  the  following  groups:  (1)  data  reduction  errors,  (2)  mean 
anomaly  errors,  (3)  representation  errors,  and  (4)  errors  of  data  reduction 
effects  on  the  estimated  quantities. 

In  Sunkel  (1981),  the  accuracies  of  the  disturbance  vector  components  at 
high  altitude  are  estimated  from  a  given  set  of  free-air  mean  anomalies  and 
their  accuracies  at  sea  level.  The  general  conclusion  from  the  study  is  that 
the  radial  component  can  be  estimated  with  *1  mgal  accuracy  at  50,000  ft 
altitude  from  the  given,  reasonable,  data  set.  To  obtain  the  same  accuracy  at 
30,000  ft,  the  accuracy  of  the  data  (especially  the  5’  *  5’  mean  anomaly  field) 
has  to  be  increased  by  60%.  For  the  horizontal  components,  the  best  given 
data  distribution  yields  only  *2.3  mgal  accuracy  at  30,000  ft  (O'.  5  in  the 
direction  of  the  gravity  vector).  for  the  achievement  of  *1  mgal,  the  given 
block  sizes  must  be  reduced  by  a  factor  of  2  out  to  a  spherical  distance  of 
30*,  and  the  overall  data  errors  reduced  by  30%.  Two  methods,  the  Poisson 
upward  continuation  integral  and  least  squares  collocation,  were  used.  The 
results  agree  within  10%. 

The  upward  continuation  of  gravity  anomalies  was  recently  investigated  by 
Cruz  and  Laakowski  (1984).  Three  procedures  were  used  and  the  results 
intercompared:  (1)  the  Poisson  integral  using  uncorrected  anomalies,  (2)  the 

Poisson  integration  using  terrain  corrected  surface  anomalies,  and  (3)  a 
procedure  called  "indirect  method"  where  the  given  anomaly  field  is  divided 
into  three  frequency  ranges.  The  three  ranges  (low,  medium,  and  high)  are 
treated  separately.  Values  at  altitudes  of  30,  10,  and  5  km  were  computed 
along  test  profiles.  The  profiles  obtained  from  terrain-uncorrected  anomalies 
have  a  negative  bias  of  0.6,  0.5,  and  0.7  mgal,  respective  to  the  above 
altitudes,  versus  the  results  from  terrain  corrected  data.  The  stTndard 
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deviation  of  the  differences  among  the  three  methods  are  on  the  order  of  0.5, 
0.6,  and  1.3  mgal  at  30,  10,  and  5  km  elevations.  A  FFT  was  also  tested  and 
the  results  agree  with  Poisson  integration  on  the  level  of  0.1  and  0.3  mgal  at 
30  and  10  km  elevations. 

With  reference  to  the  various  reductions  of  gravity  data  and  to  the 
studies  regarding  the  gravity  disturbance  vector  at  high  altitudes,  it  is 
recommended  to  perform  the  following  study: 

(1)  Considering  that  digital  terrain  models  (DTMs)  are  already  available  in  a 
gridded  form  for  many  regional  and  local  areas,  and  harmonic  coefficients 
of  the  potential  up  to  order  and  degree  180,  and  250  are  also  available, 
the  feasibility  of  obtaining  free-air  mean  anomalies  for  small  blocks  (5’  * 
5',  15*  >  15',  30’  x  30')  from  the  above  information  should  be  investigated. 

(2)  The  various  steps  of  data  reductions  and  the  feasible  methods  for 
prediction  of  the  gravity  disturbance  vector  at  high  altitudes  30,000  to 
200,000  ft)  should  be  evaluated  by  numerical  analyses.  For  the  optimal 
and  feasible  procedure,  an  algorithm  Bhould  be  prepared  consisting  of  a 
series  of  existing  or  new  subroutines  for  the  various  steps  of  the 
procedures.  The  computations  should  provide  the  accuracies  of  each  step 
involved  and  the  errors  of  the  final  values,  i.e.  of  the  vector  components 
of  the  gravity  disturbance  at  the  altitude. 

6.5.  Airborne  Gradiometry. 

In  section  5.2.6.,  two  recent  studies  by  Jekeli  (1983,  1984)  and  two  earlier 
reports  by  K.P.  Schwarz  (1976,1977)  on  this  subject  are  described.  Each 
author  prepared  a  study  on  the  accuracies  achievable  with  an  assumed 
instrumentation,  measuring  accuracies,  and  observation  pattern.  The  second 
analyses  by  the  authors  are  concerned  with  techniques  for  processing  of  the 
large  number  of  measurements. 

Jekeli  (1983)  evaluates  the  accuracy  of  the  point  gravity  vector  from  the 
five  independent  components  of  the  gradient  tenser  observed  at  altitude.  A 
uniform  bidirectional  track  spacing  at  600  m  altitude  was  used.  The 
gradiometer  accuracy  was  assumed  to  be  1.9E  standard  error.  The  accuracy 
of  the  obtained  gravity  disturbance  at  5  km  track  spacing  was  about  0.4  mgal, 
at  10  km  spacing  about  0.7  mgal,  and  at  15  km  spacing  about  1  mgal.  The 


deflection  components  were  accurate  to  016,  Oil,  and  O'. 22  respectively.  The 
change  of  the  accuracies  in  the  function  of  altitude  is  also  give  in  section 
5.2.6.  The  accuracies  refer  to  a  reference  field  corresponding  to  a  harmonic 
expansion  of  about  degree  80;  therefore,  the  figures  in  an  absolute  sense,  will 
be  downgraded  by  the  errors  of  such  field.  The  method  of  analysis  was  the 
least  squares  collocation. 

In  Jekeli  (1984)  the  possible  arrangements  for  the  auto-covariance 
functions  (Toeplitz  structures)  to  reduce  the  required  operations  for  the 
inversion  are  discussed.  These  arrangements  (possible  under  certain 
restrictions  for  the  pattern  of  the  observations),  will  not  affect  the  optimal 
solution.  Several  "non-optimal"  and  "virtually  optimal  solutions"  are  also 
discussed.  The  various  processing  techniques  are  intercompared  in  a 
numerical  analysis. 

The  accuracy  study  by  Schwarz  (1976)  investigates  the  achievable 
accuracies  of  mean  anomalies  for  various  block  sizes.  The  observed  quantities 
are  the  first  and  second  order  gradients  of  the  potential,  measured  by  a 
combined  accelerometer-gradiometer  system  (six  measurements  at  each  point). 
The  inertial  and  gravitational  forces  can  be  separated  from  the  output  of  the 
accelerometers  by  a  system  of  linear  second-order  differential  equations 
interrelating  the  gravity  disturbance  vector  with  the  observed  Becond  order 
gradients  and  the  accelerometer's  output.  By  the  integration  of  the  first 
order  gradients  along  the  flight  profile  gravity  anomalies  (Ag)  were  obtained. 
The  measurements  used  were  Ag  and  the  five  components  of  the  second  order 
gradients  in  a  4,X,r  system.  The  least  squares  collocation  technique  was  used 
for  the  analysis.  At  equal  distance  between  two  flight  profiles,  interpolation 
points  were  computed  using  a  range  of  observations  along  both  flight  profiles. 
From  the  interpolated  profiles,  the  mean  anomalies  were  computed.  With  0*.3 
spacing  of  east-west  profiles,  *3  mgal  accuracy  was  obtained  for  5'  *  5’  and 
15’  *  15’  blocks.  With  1*  spacing,  *5  mgal  was  the  accuracy  of  the  1*  *  1* 
blocks.  The  flight  altitude  was  10  km. 

The  simulation  study  of  Schwarz  (1977)  consisted  of  four  steps:  (1) 
gravity  and  gradiometry  data  were  generated  at  ground  and  flight  level,  (2) 
the  flight  level  data  were  corrupted  by  an  error  model,  (3)  the  ground  level 
anomalies  are  estimated  form  the  flight  level  data,  and  (4)  estimated  and 
generated  ground  level  data  were  compared.  Again  the  least  squares 
collocation  was  used  for  the  analysis.  Parallel  observation  profiles  with  a 


constant  intervals  for  the  observations  points  (2.6  km)  was  the  pattern  of  the 
measurements.  This  permitted  computation  of  a  "moving  operator"  (from  the 
cross-covariance  between  signal  and  observation  and  from  the  auto-covariance 
matrix  of  the  observations).  This  unchanged  operator  was  moved  along  the 
flight  profiles  from  one  set  of  data  points  to  the  next.  Th»  arrangement 
allowed  an  efficient  processing.  The  standard  errors  of  the  observed 
anomalies  and  gradients  were  assumed  to  be  *1  mgal  and  *1E  respectively. 
The  standard  error  of  a  point  along  the  interpolated  profile  on  the  ground 
was  3.3  mgal.  Comparing  the  mean  anomaly  accuracies  obtained  from  the 
accuracy  study  and  simulation,  the  average  of  the  differences  is  0.4  mgal. 

Recommendation:  The  conclusions  of  the  experimental  work  and  theoretical 
studies  in  the  area  of  airborne  gravimetry  conducted  by  AFGL  and  OSU  in  the 
late  1960’s  and  early  1970’s  was  that  the  desired  resolution  and  accuracy 
cannot  be  attained  using  accelerometer-type  sensors  on  a  moving  platform. 
The  reason  for  this  was  that  the  gravimeters  or  accelerometers  measured  not 
only  gravity  but  also  inertial  disturbances,  due  to  the  various  aircraft 
motions.  According  to  the  principle  of  equivalence,  the  two  forces  cannot  be 
rigorously  separated.  Attempts  to  separate  by  the  use  of  the  frequency 
differences  between  gravitation  and  inertia  was  only  partially  successful.  A 
rigorous  separation  is  possible  in  case  of  second  order  gradients  and  linear 
inertial  accelerations.  This  experience  initiated  the  research  to  produce  a 
moving  base  gravity  gradiometer.  The  Global  Positioning  System  (GPS) 
currently  under  full  development  changes  the  above  situation.  The  GPS  can 
provide  instantaneously  the  position  and  velocity  of  a  moving  object  with  high 
accuracy.  If  an  aircraft  equipped  with  accelerometers  measuring  the  combined 
gravitational  and  inertial  forces,  the  GPS  system  can  furnish  independently 
the  inertial  part,  i.e.  no  second  order  gradient  sensors  are  required.  This 
would  greatly  simplify  the  process  of  airborne  gravimetry.  It  is  recommended 
that  a  study  be  performed  for  the  determination  of  the  feasibility,  system 
parameters,  and  expected  accuracies  of  this  approach. 
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